From 7e5b78f72f969797a29faa1a3b36ff0b5664a1df Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 11 Mar 2021 16:25:30 +0100 Subject: [PATCH] Implementation of 3D cartesian coordinates completed. Last commit before testing. --- src/makefile | 1 + .../mesh/3DCart/moduleMesh3DCartRead.f90 | 157 +++++++++++++++++- src/modules/moduleInput.f90 | 5 + src/modules/moduleSolver.f90 | 62 +++++++ 4 files changed, 223 insertions(+), 2 deletions(-) diff --git a/src/makefile b/src/makefile index e3801fd..7bb277a 100644 --- a/src/makefile +++ b/src/makefile @@ -4,6 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \ + $(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o \ $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \ $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \ $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o \ diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 index a2197a4..40f81d0 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 @@ -9,10 +9,10 @@ MODULE moduleMesh3DCartRead END TYPE - INTERFACE connect + INTERFACE connected MODULE PROCEDURE connectedVolVol, connectedVolEdge - END INTERFACE connect + END INTERFACE connected CONTAINS !Init mesh @@ -388,6 +388,7 @@ MODULE moduleMesh3DCartRead END IF END IF + !Check surface 4 IF (.NOT. ASSOCIATED(elemA%e3)) THEN IF ((elemA%n1%n == elemB%n1%n .AND. & @@ -454,6 +455,134 @@ MODULE moduleMesh3DCartRead CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA 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 SUBROUTINE constructGlobalK(K, elem) @@ -465,6 +594,30 @@ MODULE moduleMesh3DCartRead INTEGER:: nNodes, i, j 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 MODULE moduleMesh3DCartRead diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 13d43e3..30fbc57 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -653,6 +653,7 @@ MODULE moduleInput !Read the geometry (mesh) for the case SUBROUTINE readGeometry(config) USE moduleMesh + USE moduleMesh3DCartRead, ONLY: mesh3DCartGeneric USE moduleMesh2DCylRead, ONLY: mesh2DCylGeneric USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric @@ -670,6 +671,10 @@ MODULE moduleInput !Selects the type of geometry. CALL config%get('geometry.type', geometryType, found) SELECT CASE(geometryType) + CASE ("3DCart") + !Creates a 3D cylindrical mesh + ALLOCATE(mesh3DCartGeneric:: mesh) + CASE ("2DCyl") !Creates a 2D cylindrical mesh ALLOCATE(mesh2DCylGeneric:: mesh) diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 4b0124e..d6b10a5 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -162,6 +162,68 @@ MODULE moduleSolver 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 PURE SUBROUTINE push2DCylNeutral(part, tauIn) USE moduleSpecies