Implementation of 3D cartesian coordinates completed. Last commit before

testing.
This commit is contained in:
Jorge Gonzalez 2021-03-11 16:25:30 +01:00
commit 7e5b78f72f
4 changed files with 223 additions and 2 deletions

View file

@ -4,6 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \
$(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o \
$(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \ $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \
$(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \ $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \
$(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o \ $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o \

View file

@ -9,10 +9,10 @@ MODULE moduleMesh3DCartRead
END TYPE END TYPE
INTERFACE connect INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connect END INTERFACE connected
CONTAINS CONTAINS
!Init mesh !Init mesh
@ -388,6 +388,7 @@ MODULE moduleMesh3DCartRead
END IF END IF
END IF END IF
!Check surface 4 !Check surface 4
IF (.NOT. ASSOCIATED(elemA%e3)) THEN IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. & IF ((elemA%n1%n == elemB%n1%n .AND. &
@ -454,6 +455,134 @@ MODULE moduleMesh3DCartRead
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB
!Check surface 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n1%n)) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) .OR. &
(elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n)) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
END IF
!Check surface 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) .OR. &
(elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n)) THEN
elemA%e2 => elemB
elemB%e2 => elemA
END IF
END IF
!Check surface 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF ((elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) .OR. &
(elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n)) THEN
elemA%e3 => elemB
elemB%e2 => elemA
END IF
END IF
!Check surface 4
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e4 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n1%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) .OR. &
(elemA%n1%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n)) THEN
elemA%e4 => elemB
elemB%e2 => elemA
END IF
END IF
END SUBROUTINE connectedTetraEdge END SUBROUTINE connectedTetraEdge
SUBROUTINE constructGlobalK(K, elem) SUBROUTINE constructGlobalK(K, elem)
@ -465,6 +594,30 @@ MODULE moduleMesh3DCartRead
INTEGER:: nNodes, i, j INTEGER:: nNodes, i, j
INTEGER, ALLOCATABLE:: n(:) INTEGER, ALLOCATABLE:: n(:)
SELECT TYPE(elem)
TYPE IS(meshVol3DCartTetra)
nNodes = 4
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
CLASS DEFAULT
nNodes = 0
ALLOCATE(localK(1:1, 1:1))
localK = 0.D0
ALLOCATE(n(1:1))
n = 0
END SELECT
DO i = 1, nNodes
DO j = 1, nNodes
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK END SUBROUTINE constructGlobalK
END MODULE moduleMesh3DCartRead END MODULE moduleMesh3DCartRead

View file

@ -653,6 +653,7 @@ MODULE moduleInput
!Read the geometry (mesh) for the case !Read the geometry (mesh) for the case
SUBROUTINE readGeometry(config) SUBROUTINE readGeometry(config)
USE moduleMesh USE moduleMesh
USE moduleMesh3DCartRead, ONLY: mesh3DCartGeneric
USE moduleMesh2DCylRead, ONLY: mesh2DCylGeneric USE moduleMesh2DCylRead, ONLY: mesh2DCylGeneric
USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric
USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric
@ -670,6 +671,10 @@ MODULE moduleInput
!Selects the type of geometry. !Selects the type of geometry.
CALL config%get('geometry.type', geometryType, found) CALL config%get('geometry.type', geometryType, found)
SELECT CASE(geometryType) SELECT CASE(geometryType)
CASE ("3DCart")
!Creates a 3D cylindrical mesh
ALLOCATE(mesh3DCartGeneric:: mesh)
CASE ("2DCyl") CASE ("2DCyl")
!Creates a 2D cylindrical mesh !Creates a 2D cylindrical mesh
ALLOCATE(mesh2DCylGeneric:: mesh) ALLOCATE(mesh2DCylGeneric:: mesh)

View file

@ -162,6 +162,68 @@ MODULE moduleSolver
END SUBROUTINE doPushes END SUBROUTINE doPushes
!Push neutral particles in 3D cartesian coordinates
PURE SUBROUTINE push3DCartNeutral(part, tauIn)
USE moduleSPecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
part_temp = part
!x
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!y
part_temp%v(2) = part%v(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
!z
part_temp%v(3) = part%v(3)
part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push3DCartNeutral
!Push charged particles in 3D cartesian coordinates
PURE SUBROUTINE push3DCartCharged(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: qmEFt(1:3)
part_temp = part
!Get the electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
!x
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!y
part_temp%v(2) = part%v(2) + qmEFt(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
!z
part_temp%v(3) = part%v(3) + qmEFt(3)
part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push3DCartCharged
!Push one particle. Boris pusher for 2D Cyl Neutral particle !Push one particle. Boris pusher for 2D Cyl Neutral particle
PURE SUBROUTINE push2DCylNeutral(part, tauIn) PURE SUBROUTINE push2DCylNeutral(part, tauIn)
USE moduleSpecies USE moduleSpecies