fpakc/src/modules/mesh/1DRad/moduleMesh1DRad.f90
JGonzalez 4e9514876e Number of collisions per collision pair
Now the number of collisions is calculated per species pair. This allows
that the randomly particles selected for collisions do not have
collisions assigned.
2022-04-23 19:00:33 +02:00

638 lines
16 KiB
Fortran

!moduleMesh1DRad: 1D radial module
! x == r
! y == theta (unused)
! z == unused
MODULE moduleMesh1DRad
USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE
REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad
!Element coordinates
REAL(8):: r = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initNode1DRad
PROCEDURE, PASS:: getCoordinates => getCoord1DRad
END TYPE meshNode1DRad
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge1DRad
!Element coordinates
REAL(8):: r = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS
PROCEDURE, PASS:: init => initEdge1DRad
PROCEDURE, PASS:: getNodes => getNodes1DRad
PROCEDURE, PASS:: intersection => intersection1DRad
PROCEDURE, PASS:: randPos => randPos1DRad
END TYPE meshEdge1DRad
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad
CONTAINS
PROCEDURE, PASS:: detJac => detJ1DRad
PROCEDURE, PASS:: invJac => invJ1DRad
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol1DRad
ABSTRACT INTERFACE
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
END FUNCTION dPsi_interface
PURE SUBROUTINE partialDer_interface(self, dPsi, dx)
IMPORT meshVol1DRad
CLASS(meshVol1DRad), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1):: dx
END SUBROUTINE partialDer_interface
END INTERFACE
TYPE, PUBLIC, EXTENDS(meshVol1DRad):: meshVol1DRadSegm
!Element coordinates
REAL(8):: r(1:2)
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
!Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS
PROCEDURE, PASS:: init => initVol1DRadSegm
PROCEDURE, PASS:: randPos => randPos1DRadSeg
PROCEDURE, PASS:: area => areaRad
PROCEDURE, NOPASS:: fPsi => fPsiRad
PROCEDURE, NOPASS:: dPsi => dPsiRad
PROCEDURE, PASS:: partialDer => partialDerRad
PROCEDURE, PASS:: elemK => elemKRad
PROCEDURE, PASS:: elemF => elemFRad
PROCEDURE, NOPASS:: inside => insideRad
PROCEDURE, PASS:: gatherEF => gatherEFRad
PROCEDURE, PASS:: gatherMF => gatherMFRad
PROCEDURE, PASS:: getNodes => getNodesRad
PROCEDURE, PASS:: phy2log => phy2logRad
PROCEDURE, PASS:: nextElement => nextElementRad
END TYPE meshVol1DRadSegm
CONTAINS
!NODE FUNCTIONS
!Init node element
SUBROUTINE initNode1DRad(self, n, r)
USE moduleSpecies
USE moduleRefParam
USE OMP_LIB
IMPLICIT NONE
CLASS(meshNode1DRad), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%r = r(1)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output
ALLOCATE(self%output(1:nSpecies))
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initNode1DRad
PURE FUNCTION getCoord1DRad(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNode1DRad), INTENT(in):: self
REAL(8):: r(1:3)
r = (/ self%r, 0.D0, 0.D0 /)
END FUNCTION getCoord1DRad
!EDGE FUNCTIONS
!Inits edge element
SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface)
USE moduleSpecies
USE moduleBoundary
USE moduleErrors
IMPLICIT NONE
CLASS(meshEdge1DRad), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
REAL(8), DIMENSION(1:3):: r1
INTEGER:: s
self%n = n
self%n1 => mesh%nodes(p(1))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
self%r = r1(1)
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
self%normal = self%normal/NORM2(self%normal)
!Boundary index
self%boundary => boundary(bt)
ALLOCATE(self%fboundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
CALL pointBoundaryFunction(self, s)
END DO
!Physical Surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdge1DRad
!Get nodes from edge
PURE FUNCTION getNodes1DRad(self) RESULT(n)
IMPLICIT NONE
CLASS(meshEdge1DRad), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1))
n = (/ self%n1%n /)
END FUNCTION getNodes1DRad
PURE FUNCTION intersection1DRad(self, r0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge1DRad), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0
REAL(8), DIMENSION(1:3):: r
r = (/ self%r, 0.D0, 0.D0 /)
END FUNCTION intersection1DRad
!Calculates a 'random' position in edge
FUNCTION randPos1DRad(self) RESULT(r)
CLASS(meshEdge1DRad), INTENT(in):: self
REAL(8):: r(1:3)
r = (/ self%r, 0.D0, 0.D0 /)
END FUNCTION randPos1DRad
!VOLUME FUNCTIONS
!SEGMENT FUNCTIONS
!Init segment element
SUBROUTINE initVol1DRadSegm(self, n, p, nodes)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
REAL(8), DIMENSION(1:3):: r1, r2
self%n = n
self%n1 => nodes(p(1))%obj
self%n2 => nodes(p(2))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
self%r = (/ r1(1), r2(1) /)
!Assign node volume
CALL self%area()
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies))
ALLOCATE(self%totalWeight(1:nSpecies))
END SUBROUTINE initVol1DRadSegm
!Calculates a random position in 1D volume
FUNCTION randPos1DRadSeg(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
xii(1) = random(-1.D0, 1.D0)
xii(2:3) = 0.D0
fPsi = self%fPsi(xii)
r(1) = DOT_PRODUCT(fPsi, self%r)
END FUNCTION randPos1DRadSeg
!Computes element area
PURE SUBROUTINE areaRad(self)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(inout):: self
REAL(8):: l !element length
REAL(8):: fPsi(1:2)
REAL(8):: r
REAL(8):: detJ
REAL(8):: Xii(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!1 point Gauss integral
Xii = 0.D0
fPsi = self%fPsi(Xii)
detJ = self%detJac(Xii)
r = DOT_PRODUCT(fPsi, self%r)
l = 2.D0*detJ
self%volume = r*l
self%arNodes = fPsi*r*l
END SUBROUTINE areaRad
!Computes element functions at point xii
PURE FUNCTION fPsiRad(xi) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:2))
fPsi(1) = 1.D0 - xi(1)
fPsi(2) = 1.D0 + xi(1)
fPsi = fPsi * 5.D-1
END FUNCTION fPsiRad
!Computes element derivative shape function at Xii
PURE FUNCTION dPsiRad(xi) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:1, 1:2))
dPsi(1, 1) = -5.D-1
dPsi(1, 2) = 5.D-1
END FUNCTION dPsiRad
!Computes partial derivatives of coordinates
PURE SUBROUTINE partialDerRad(self, dPsi, dx)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1):: dx
dx(1) = DOT_PRODUCT(dPsi(1,:), self%r)
END SUBROUTINE partialDerRad
!Computes local stiffness matrix
PURE FUNCTION elemKRad(self) RESULT(localK)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), ALLOCATABLE:: localK(:,:)
REAL(8):: Xii(1:3)
REAL(8):: dPsi(1:1, 1:2)
REAL(8):: invJ(1), detJ
REAL(8):: r, fPsi(1:2)
INTEGER:: l
ALLOCATE(localK(1:2, 1:2))
localK = 0.D0
Xii = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
dPsi = self%dPsi(Xii)
detJ = self%detJac(Xii, dPsi)
invJ = self%invJac(Xii, dPsi)
fPsi = self%fPsi(Xii)
r = DOT_PRODUCT(fPsi, self%r)
localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
r*wSeg(l)/detJ
END DO
localK = localK*PI2
END FUNCTION elemKRad
PURE FUNCTION elemFRad(self, source) RESULT(localF)
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:2)
REAL(8):: detJ, f, r
REAL(8):: Xii(1:3)
INTEGER:: l
ALLOCATE(localF(1:2))
localF = 0.D0
Xii = 0.D0
DO l = 1, 3
xii(1) = corSeg(l)
detJ = self%detJac(Xii)
fPsi = self%fPsi(Xii)
r = DOT_PRODUCT(fPsi, self%r)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*r*wSeg(l)*detJ
END DO
END FUNCTION elemFRad
PURE FUNCTION insideRad(xi) RESULT(ins)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
LOGICAL:: ins
ins = xi(1) >=-1.D0 .AND. &
xi(1) <= 1.D0
END FUNCTION insideRad
!Gathers EF at position Xii
PURE FUNCTION gatherEFRad(self, xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1, 1:2)
REAL(8):: phi(1:2)
REAL(8):: EF(1:3)
REAL(8):: invJ
phi = (/ self%n1%emData%phi, &
self%n2%emData%phi /)
dPsi = self%dPsi(xi)
invJ = self%invJac(xi, dPsi)
EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ
EF(2) = 0.D0
EF(3) = 0.D0
END FUNCTION gatherEFRad
PURE FUNCTION gatherMFRad(self, xi) RESULT(MF)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: fPsi(1:2)
REAL(8):: MF_Nodes(1:2, 1:3)
REAL(8):: MF(1:3)
REAL(8):: invJ
MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), &
self%n2%emData%B(1) /)
MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), &
self%n2%emData%B(2) /)
MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), &
self%n2%emData%B(3) /)
fPsi = self%fPsi(xi)
MF = MATMUL(fPsi, MF_Nodes)
END FUNCTION gatherMFRad
!Get nodes from 1D volume
PURE FUNCTION getNodesRad(self) RESULT(n)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:2))
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesRad
PURE FUNCTION phy2logRad(self, r) RESULT(xN)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
xN = 0.D0
xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0
END FUNCTION phy2logRad
!Get next element for a logical position xi
SUBROUTINE nextElementRad(self, xi, nextElement)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(meshElement), POINTER, INTENT(out):: nextElement
NULLIFY(nextElement)
IF (xi(1) < -1.D0) THEN
nextElement => self%e2
ELSEIF (xi(1) > 1.D0) THEN
nextElement => self%e1
END IF
END SUBROUTINE nextElementRad
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes the element Jacobian determinant
PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVol1DRad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dJ
REAL(8):: dx(1)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
dJ = dx(1)
END FUNCTION detJ1DRad
!Computes the invers Jacobian
PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVol1DRad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1)
REAL(8):: invJ
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
invJ = 1.D0/dx(1)
END FUNCTION invJ1DRad
SUBROUTINE connectMesh1DRad(self)
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout):: self
INTEGER:: e, et
DO e = 1, self%numVols
!Connect Vol-Vol
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
SELECT TYPE(self)
TYPE IS(meshParticles)
!Connect Vol-Edge
DO et = 1, self%numEdges
CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
END DO
END SELECT
END DO
END SUBROUTINE connectMesh1DRad
SUBROUTINE connectVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVol1DRadSegm)
SELECT TYPE(elemB)
TYPE IS(meshVol1DRadSegm)
CALL connectSegmSegm(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectVolVol
SUBROUTINE connectSegmSegm(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END SUBROUTINE connectSegmSegm
SUBROUTINE connectVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS (meshVol1DRadSegm)
SELECT TYPE(elemB)
CLASS IS(meshEdge1DRad)
CALL connectSegmEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectVolEdge
SUBROUTINE connectSegmEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol1DRadSegm), INTENT(inout), TARGET:: elemA
CLASS(meshEdge1DRad), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END SUBROUTINE connectSegmEdge
END MODULE moduleMesh1DRad