From 7e5b78f72f969797a29faa1a3b36ff0b5664a1df Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 11 Mar 2021 16:25:30 +0100 Subject: [PATCH 1/4] 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 From c236c5e0e26dcd4ab3046527969ca84ddaf7bf00 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 15 Mar 2021 10:00:34 +0100 Subject: [PATCH 2/4] Fixing some issues in 3D Cartesian coordinates. Included 3D pusher. Still there are issues linking a volume to an edge. --- .../mesh/1DCart/moduleMesh1DCartRead.f90 | 3 + .../mesh/1DRad/moduleMesh1DRadRead.f90 | 3 + .../mesh/2DCart/moduleMesh2DCartRead.f90 | 21 ++ .../mesh/2DCyl/moduleMesh2DCylRead.f90 | 18 ++ src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 2 +- .../mesh/3DCart/moduleMesh3DCartRead.f90 | 208 ++++++------------ src/modules/moduleInject.f90 | 9 +- src/modules/moduleSolver.f90 | 7 + 8 files changed, 123 insertions(+), 148 deletions(-) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 b/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 index 720e61d..d69cfdc 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCartRead.f90 @@ -217,6 +217,9 @@ MODULE moduleMesh1DCartRead elemA%e1 => elemB elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 b/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 index 1a9f25c..4f14b5f 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRadRead.f90 @@ -218,6 +218,9 @@ MODULE moduleMesh1DRadRead 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. & diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 index b91627a..b0f2ee3 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 @@ -439,6 +439,9 @@ MODULE moduleMesh2DCartRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -455,6 +458,9 @@ MODULE moduleMesh2DCartRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -471,6 +477,9 @@ MODULE moduleMesh2DCartRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB @@ -487,6 +496,9 @@ MODULE moduleMesh2DCartRead elemA%e4 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n4%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e4 => elemB @@ -511,6 +523,9 @@ MODULE moduleMesh2DCartRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -527,6 +542,9 @@ MODULE moduleMesh2DCartRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -543,6 +561,9 @@ MODULE moduleMesh2DCartRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 index 085c3c3..379a0fb 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 @@ -439,6 +439,9 @@ MODULE moduleMesh2DCylRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -455,6 +458,9 @@ MODULE moduleMesh2DCylRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -471,6 +477,9 @@ MODULE moduleMesh2DCylRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB @@ -511,6 +520,9 @@ MODULE moduleMesh2DCylRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB @@ -527,6 +539,9 @@ MODULE moduleMesh2DCylRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB @@ -543,6 +558,9 @@ MODULE moduleMesh2DCylRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 1c713c9..0464c7f 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -146,7 +146,7 @@ MODULE moduleMesh3DCart self%n = n self%n1 => mesh%nodes(p(1))%obj - self%n3 => mesh%nodes(p(2))%obj + self%n2 => mesh%nodes(p(2))%obj self%n3 => mesh%nodes(p(3))%obj !Get element coordinates r1 = self%n1%getCoordinates() diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 index 40f81d0..c6b021a 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 @@ -207,6 +207,28 @@ MODULE moduleMesh3DCartRead END SUBROUTINE connectedVolEdge + PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) + IMPLICIT NONE + + INTEGER, DIMENSION(1:3), INTENT(in):: nodesA, nodesB + LOGICAL:: coincident + INTEGER:: i + + coincident = .FALSE. + DO i = 1, 3 + IF (ANY(nodesA(i) == nodesB)) THEN + coincident = .TRUE. + + ELSE + coincident = .FALSE. + EXIT + + END IF + + END DO + + END FUNCTION coincidentNodes + SUBROUTINE connectedTetraTetra(elemA, elemB) IMPLICIT NONE @@ -217,54 +239,26 @@ MODULE moduleMesh3DCartRead !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 + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e1 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%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%n4%n .AND. & - elemA%n3%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e1 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%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%n4%n .AND. & - elemA%n3%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e1 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n3%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e1 => elemB elemB%e4 => elemA @@ -275,54 +269,26 @@ MODULE moduleMesh3DCartRead !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 + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e2 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%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%n4%n .AND. & - elemA%n4%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%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%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e4 => elemA @@ -333,54 +299,26 @@ MODULE moduleMesh3DCartRead !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 + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e3 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n2%n == elemB%n4%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%n4%n .AND. & - elemA%n4%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n2%n == elemB%n4%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%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n2%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e4 => elemA @@ -391,54 +329,26 @@ MODULE moduleMesh3DCartRead !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 + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e4 => elemB elemB%e1 => elemA - ELSEIF ((elemA%n1%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%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%n4%n .AND. & - elemA%n4%n == elemB%n2%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e4 => elemB elemB%e2 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%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%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e4 => elemB elemB%e3 => elemA - ELSEIF ((elemA%n1%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n4%n) .OR. & - (elemA%n1%n == elemB%n4%n .AND. & - elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) .OR. & - (elemA%n1%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n4%n .AND. & - elemA%n4%n == elemB%n1%n)) THEN + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e4 => elemB elemB%e4 => elemA @@ -483,8 +393,12 @@ MODULE moduleMesh3DCartRead elemA%e1 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + END IF + END IF !Check surface 2 @@ -515,6 +429,9 @@ MODULE moduleMesh3DCartRead elemA%e2 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + END IF END IF @@ -547,6 +464,9 @@ MODULE moduleMesh3DCartRead elemA%e3 => elemB elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + END IF END IF diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9a02cb1..e614cdf 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -230,12 +230,15 @@ MODULE moduleInject !Random position in edge partInj(n)%r = randomEdge%randPos() !Volume associated to the edge: - IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN + IF (ASSOCIATED(randomEdge%e1)) THEN partInj(n)%vol = randomEdge%e1%n - - ELSE + + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN partInj(n)%vol = randomEdge%e2%n + ELSE + PRINT *, "ERROR NO VOL ASSOCIATED TO EDGE" + END IF partInj(n)%v = (/ self%v(1)%obj%randomVel(), & diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index d6b10a5..c76b6c9 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -67,6 +67,13 @@ MODULE moduleSolver REAL(8):: tau, tauSp SELECT CASE(pusherType) + !3D Cartesian + CASE('3DCartNeutral') + self%pushParticle => push3DCartNeutral + + CASE('3DCartCharged') + self%pushParticle => push3DCartCharged + !2D Cylindrical CASE('2DCylNeutral') self%pushParticle => push2DCylNeutral From 8cf50cda68c1cbd49cf19c597d4b27ba975b6b7c Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 20 Mar 2021 13:08:01 +0100 Subject: [PATCH 3/4] First implementation of ionization boundary. Still some work to do. --- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 62 +++++----- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 18 +-- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 18 +-- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 25 +--- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 18 +-- src/modules/mesh/moduleMeshBoundary.f90 | 113 +++++++++++++++++++ src/modules/moduleBoundary.f90 | 45 +++++++- src/modules/moduleInject.f90 | 4 + src/modules/moduleInput.f90 | 26 +++++ src/modules/moduleList.f90 | 2 + src/modules/moduleSolver.f90 | 21 +++- 11 files changed, 247 insertions(+), 105 deletions(-) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index d45c7a0..366fa84 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -7,6 +7,9 @@ MODULE moduleMesh1DCart 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):: meshNode1DCart !Element coordinates REAL(8):: x = 0.D0 @@ -152,23 +155,7 @@ MODULE moduleMesh1DCart ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies - SELECT TYPE(obj => self%boundary%bTypes(s)%obj) - TYPE IS(boundaryAbsorption) - self%fBoundary(s)%apply => absorption - - TYPE IS(boundaryReflection) - self%fBoundary(s)%apply => reflection - - TYPE IS(boundaryTransparent) - self%fBoundary(s)%apply => transparent - - TYPE IS(boundaryWallTemperature) - self%fBoundary(s)%apply => wallTemperature - - CLASS DEFAULT - CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DCart') - - END SELECT + CALL pointBoundaryFunction(self, s) END DO @@ -322,22 +309,29 @@ MODULE moduleMesh1DCart END SUBROUTINE partialDerSegm !Computes local stiffness matrix - PURE FUNCTION elemKSegm(self) RESULT(ke) + FUNCTION elemKSegm(self) RESULT(ke) IMPLICIT NONE CLASS(meshVol1DCartSegm), INTENT(in):: self REAL(8):: ke(1:2,1:2) REAL(8):: Xii(1:3) REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ + REAL(8):: invJ(1), detJ + INTEGER:: l ke = 0.D0 - Xii = (/ 0.D0, 0.D0, 0.D0 /) - dPsi = self%dPsi(Xii) - invJ = self%invJac(Xii, dPsi) - ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) - ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) - ke = 2.D0*ke*invJ + 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) + ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + wSeg(l)/detJ + + END DO END FUNCTION elemKSegm @@ -348,14 +342,22 @@ MODULE moduleMesh1DCart REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:2) - REAL(8):: detJ + REAL(8):: detJ, f REAL(8):: Xii(1:3) + INTEGER:: l - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) ALLOCATE(localF(1:2)) - localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ + localF = 0.D0 + Xii = 0.D0 + + DO l = 1, 3 + xii(1) = corSeg(l) + detJ = self%detJac(Xii) + fPsi = self%fPsi(Xii) + f = DOT_PRODUCT(fPsi, source) + localF = localF + f*fPsi*wSeg(l)*detJ + + END DO END FUNCTION elemFSegm diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 1d2bdba..d0f8311 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -153,23 +153,7 @@ MODULE moduleMesh1DRad ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies - SELECT TYPE(obj => self%boundary%bTypes(s)%obj) - TYPE IS(boundaryAbsorption) - self%fBoundary(s)%apply => absorption - - TYPE IS(boundaryReflection) - self%fBoundary(s)%apply => reflection - - TYPE IS(boundaryTransparent) - self%fBoundary(s)%apply => transparent - - TYPE IS(boundaryWallTemperature) - self%fBoundary(s)%apply => wallTemperature - - CLASS DEFAULT - CALL criticalError("Boundary type not defined in this geometry", 'initEdge1DRad') - - END SELECT + CALL pointBoundaryFunction(self, s) END DO diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index ac22b28..40c218c 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -200,23 +200,7 @@ MODULE moduleMesh2DCart ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies - SELECT TYPE(obj => self%boundary%bTypes(s)%obj) - TYPE IS(boundaryAbsorption) - self%fBoundary(s)%apply => absorption - - TYPE IS(boundaryReflection) - self%fBoundary(s)%apply => reflection - - TYPE IS(boundaryTransparent) - self%fBoundary(s)%apply => transparent - - TYPE IS(boundaryWallTemperature) - self%fBoundary(s)%apply => wallTemperature - - CLASS DEFAULT - CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCart') - - END SELECT + CALL pointBoundaryFunction(self, s) END DO diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 8ef5c6d..e0f3a27 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -200,26 +200,7 @@ MODULE moduleMesh2DCyl ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies - SELECT TYPE(obj => self%boundary%bTypes(s)%obj) - TYPE IS(boundaryAbsorption) - self%fBoundary(s)%apply => absorption - - TYPE IS(boundaryReflection) - self%fBoundary(s)%apply => reflection - - TYPE IS(boundaryTransparent) - self%fBoundary(s)%apply => transparent - - TYPE IS(boundaryWallTemperature) - self%fBoundary(s)%apply => wallTemperature - - TYPE IS(boundaryAxis) - self%fBoundary(s)%apply => symmetryAxis - - CLASS DEFAULT - CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCyl') - - END SELECT + CALL pointBoundaryFunction(self, s) END DO @@ -459,7 +440,9 @@ MODULE moduleMesh2DCyl detJ = self%detJac(xi,dPsi) invJ = self%invJac(xi,dPsi) r = DOT_PRODUCT(fPsi,self%r) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wQuad(l)*wQuad(m)/detJ + ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + r*wQuad(l)*wQuad(m)/detJ END DO END DO diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 1c713c9..e1b8907 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -166,23 +166,7 @@ MODULE moduleMesh3DCart ALLOCATE(self%fBoundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies - SELECT TYPE(obj => self%boundary%bTypes(s)%obj) - TYPE IS(boundaryAbsorption) - self%fBoundary(s)%apply => absorption - - TYPE IS(boundaryReflection) - self%fBoundary(s)%apply => reflection - - TYPE IS(boundaryTransparent) - self%fBoundary(s)%apply => transparent - - TYPE IS(boundaryWallTemperature) - self%fBoundary(s)%apply => wallTemperature - - CLASS DEFAULT - CALL criticalError("Boundary type not defined in this geometry", 'initEdge3DCart') - - END SELECT + CALL pointBoundaryFunction(self, s) END DO diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 97602f4..2252b84 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -108,6 +108,85 @@ MODULE moduleMeshBoundary END SUBROUTINE wallTemperature + !Ionization surface: an electron will pass through the surface + ! and create an ion-electron pair based on a neutral background + SUBROUTINE ionization(edge, part) + USE moduleList + USE moduleSpecies + USE moduleMesh + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + CLASS(particle), INTENT(inout):: part + REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass + REAL(8):: ionizationRate + INTEGER:: ionizationPair, p + TYPE(particle), POINTER:: newElectron + TYPE(particle), POINTER:: newIon + + SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) + TYPE IS(boundaryIonization) + mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m) + vRel = SUM(DABS(part%v-bound%v0)) + eRel = mRel*vRel**2*5.D-1 + + IF (eRel > bound%eThreshold) THEN + !TODO: Check units + ionizationRate = bound%n0*bound%crossSection%get(eRel) + + ionizationPair = INT(ionizationRate*tauMin*ti_ref/species(bound%sp)%obj%weight) + + !Create the new pair of particles + DO p = 1, ionizationPair + ALLOCATE(newElectron) + ALLOCATE(newIon) + + newElectron%sp = part%sp + newIon%sp = bound%sp + + newElectron%v = 10.D0*bound%v0 + (1.D0 + bound%deltaV/NORM2(bound%v0)) + newIon%v = bound%v0 + + newElectron%r = edge%randPos() + newIon%r = newElectron%r + + newElectron%vol = part%vol + newIon%vol = part%vol + + newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r) + newIon%xi = newElectron%xi + + newElectron%qm = part%qm + SELECT TYPE(spe => species(bound%sp)%obj) + TYPE IS(speciesCharged) + newIon%qm = spe%qm + + END SELECT + + newElectron%weight = species(bound%sp)%obj%weight + newIon%weight = newElectron%weight + + newElectron%n_in = .TRUE. + newIon%n_in = .TRUE. + + !Add particles to list + CALL OMP_SET_LOCK(lockSurfaces) + CALL partSurfaces%add(newElectron) + CALL partSurfaces%add(newIon) + CALL OMP_UNSET_LOCK(lockSurfaces) + + END DO + + !Removes ionizing electron + part%n_in = .FALSE. + + END IF + + END SELECT + + END SUBROUTINE ionization + !Symmetry axis. Dummy function SUBROUTINE symmetryAxis(edge, part) USE moduleSpecies @@ -118,4 +197,38 @@ MODULE moduleMeshBoundary END SUBROUTINE symmetryAxis + !Points the boundary function to specific type + SUBROUTINE pointBoundaryFunction(edge, s) + USE moduleErrors + IMPLICIT NONE + + CLASS(meshEdge), INTENT(inout):: edge + INTEGER, INTENT(in):: s !Species index + + SELECT TYPE(obj => edge%boundary%bTypes(s)%obj) + TYPE IS(boundaryAbsorption) + edge%fBoundary(s)%apply => absorption + + TYPE IS(boundaryReflection) + edge%fBoundary(s)%apply => reflection + + TYPE IS(boundaryTransparent) + edge%fBoundary(s)%apply => transparent + + TYPE IS(boundaryWallTemperature) + edge%fBoundary(s)%apply => wallTemperature + + TYPE IS(boundaryIonization) + edge%fBoundary(s)%apply => ionization + + TYPE IS(boundaryAxis) + edge%fBoundary(s)%apply => symmetryAxis + + CLASS DEFAULT + CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction') + + END SELECT + + END SUBROUTINE pointBoundaryFunction + END MODULE moduleMeshBoundary diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 index 0f54ee9..66a98df 100644 --- a/src/modules/moduleBoundary.f90 +++ b/src/modules/moduleBoundary.f90 @@ -1,4 +1,5 @@ MODULE moduleBoundary + USE moduleTable !Generic type for boundaries TYPE, PUBLIC:: boundaryGeneric @@ -24,7 +25,7 @@ MODULE moduleBoundary END TYPE boundaryTransparent - !Transparent boundary + !Wall Temperature boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature !Thermal velocity of the wall: square root(Wall temperature X specific heat) REAL(8):: vTh @@ -32,6 +33,18 @@ MODULE moduleBoundary END TYPE boundaryWallTemperature + !Ionization boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization + REAL(8):: m0, n0, v0(1:3), T0 !Properties of background neutrals. + INTEGER:: sp !Ion species + TYPE(table1D):: crossSection + REAL(8):: eThreshold + REAL(8):: deltaV + + CONTAINS + + END TYPE boundaryIonization + !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis CONTAINS @@ -86,4 +99,34 @@ MODULE moduleBoundary END SUBROUTINE initWallTemperature + SUBROUTINE initIonization(boundary, m0, n0, v0, T0, speciesID, crossSection, eThreshold) + USE moduleRefParam + USE moduleSpecies + USE moduleCaseParam + USE moduleConstParam + IMPLICIT NONE + + CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary + REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 + INTEGER:: speciesID + CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection + REAL(8), INTENT(in):: eThreshold + + ALLOCATE(boundaryIonization:: boundary) + + SELECT TYPE(boundary) + TYPE IS(boundaryIonization) + boundary%m0 = m0 / m_ref + boundary%n0 = n0 + boundary%v0 = v0 / v_ref + boundary%T0 = T0 / T_ref + boundary%sp = speciesID + CALL boundary%crossSection%init(crossSection) + CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) + boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2) + + END SELECT + + END SUBROUTINE initIonization + END MODULE moduleBoundary diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9a02cb1..becf0a6 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -98,6 +98,10 @@ MODULE moduleInject !Input current in Ampers self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight)) + CASE ("part/s") + !Input current in Ampers + self%nParticles = INT(flow*tauMin*ti_ref/species(sp)%obj%weight) + CASE DEFAULT CALL criticalError("No support for units: " // units, 'initInject') diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 13d43e3..897ddb2 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -599,6 +599,13 @@ MODULE moduleInput CHARACTER(2):: istring, sString CHARACTER(:), ALLOCATABLE:: object, bType REAL(8):: Tw, cw !Wall temperature and specific heat + !Neutral Properties + REAL(8):: m0, n0 + REAL(8), DIMENSION(:), ALLOCATABLE:: v0 + REAL(8):: T0 + REAL(8):: eThreshold !Energy threshold + INTEGER:: speciesID + CHARACTER(:), ALLOCATABLE:: speciesName, crossSection LOGICAL:: found INTEGER:: nTypes @@ -628,6 +635,25 @@ MODULE moduleInput CASE('transparent') ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj) + CASE('ionization') + CALL config%get(object // '.neutral.name', speciesName, found) + IF (.NOT. found) CALL criticalError("missing parameter 'name' for neutrals in ionization", 'readBoundary') + speciesID = speciesName2Index(speciesName) + CALL config%get(object // '.neutral.mass', m0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'mass' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.neutral.density', n0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.neutral.velocity', v0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.neutral.temperature', T0, found) + IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary') + CALL config%get(object // '.energyThreshold', eThreshold, found) + IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary') + CALL config%get(object // '.crossSection', crossSection, found) + IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary') + + CALL initIonization(boundary(i)%bTypes(s)%obj, m0, n0, v0, T0, speciesID, crossSection, eThreshold) + CASE('wallTemperature') CALL config%get(object // '.temperature', Tw, found) IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary') diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index 15c4b23..c08061a 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -24,6 +24,8 @@ MODULE moduleList INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles TYPE(listNode):: partCollisions !Particles created in collisional process INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles + TYPE(listNode):: partSurfaces !Particles created in surface interactions + INTEGER(KIND=OMP_LOCK_KIND):: lockSurfaces !Lock for the NA list of particles TYPE(listNode):: partInitial !Initial distribution of particles TYPE pointerArray diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 4b0124e..d75f073 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -437,7 +437,7 @@ MODULE moduleSolver INTEGER:: nn, n, e INTEGER, SAVE:: nPartNew - INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions + INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) TYPE(lNode), POINTER:: partCurr, partNext @@ -458,13 +458,15 @@ MODULE moduleSolver nWScheme = partWScheme%amount !$OMP SECTION nCollisions = partCollisions%amount + !$OMP SECTION + nSurfaces = partSurfaces%amount !$OMP END SECTIONS !$OMP BARRIER !$OMP SINGLE CALL MOVE_ALLOC(partOld, partTemp) - nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces ALLOCATE(partOld(1:nPartNew)) !$OMP END SINGLE @@ -522,6 +524,21 @@ MODULE moduleSolver IF (ASSOCIATED(partCollisions%tail)) NULLIFY(partCollisions%tail) partCollisions%amount = 0 + !$OMP SECTION + !Reset particles from surface process + nn = nInjIn + nOldIn + nWScheme + nCollisions + partCurr => partSurfaces%head + DO n = 1, nSurfaces + partNext => partCurr%next + partOld(nn+n) = partCurr%part + DEALLOCATE(partCurr) + partCurr => partNext + + END DO + IF (ASSOCIATED(partSurfaces%head)) NULLIFY(partSurfaces%head) + IF (ASSOCIATED(partSurfaces%tail)) NULLIFY(partSurfaces%tail) + partSurfaces%amount = 0 + !$OMP SECTION !Reset output in nodes DO e = 1, mesh%numNodes From db6b0a2c03128352af55a0abf9b102304a57125d Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 22 Mar 2021 12:39:34 +0100 Subject: [PATCH 4/4] Fixed an issue with reflection of particles in all geometries and also assigning the normal vector in 2D and 3D. 3D Cartesian geometry is working properly, although it needs testing. Still issue with ionization boundary. --- src/makefile | 2 +- src/modules/makefile | 6 +- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 7 +- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 56 +++-- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 25 ++- .../mesh/2DCart/moduleMesh2DCartRead.f90 | 70 +++--- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 26 ++- .../mesh/2DCyl/moduleMesh2DCylRead.f90 | 90 +++++--- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 41 ++-- .../mesh/3DCart/moduleMesh3DCartRead.f90 | 211 +++++++++--------- src/modules/mesh/moduleMesh.f90 | 4 +- src/modules/mesh/moduleMeshBoundary.f90 | 20 +- src/modules/moduleInject.f90 | 3 +- src/modules/moduleMath.f90 | 42 ++++ src/modules/moduleOutput.f90 | 11 +- 15 files changed, 349 insertions(+), 265 deletions(-) create mode 100644 src/modules/moduleMath.f90 diff --git a/src/makefile b/src/makefile index 7bb277a..b398134 100644 --- a/src/makefile +++ b/src/makefile @@ -3,7 +3,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ $(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)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ $(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o \ $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \ $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \ diff --git a/src/modules/makefile b/src/modules/makefile index 28fc4ec..cea6d73 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -2,11 +2,11 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \ moduleOutput.o moduleInput.o moduleSolver.o \ moduleCollisions.o moduleTable.o moduleParallel.o \ - moduleEM.o moduleRandom.o + moduleEM.o moduleRandom.o moduleMath.o all: $(OBJS) mesh.o -mesh.o: moduleCollisions.o moduleBoundary.o +mesh.o: moduleCollisions.o moduleBoundary.o moduleMath.o $(MAKE) -C mesh all moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90 @@ -21,7 +21,7 @@ moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90 moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90 +moduleOutput.o: moduleMath.o moduleSpecies.o moduleRefParam.o moduleOutput.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ moduleRandom.o: moduleConstParam.o moduleRandom.f90 diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 366fa84..020b971 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -148,7 +148,6 @@ MODULE moduleMesh1DCart self%x = r1(1) self%normal = (/ 1.D0, 0.D0, 0.D0 /) - self%normal = self%normal/NORM2(self%normal) !Boundary index self%boundary => boundary(bt) @@ -176,11 +175,11 @@ MODULE moduleMesh1DCart END FUNCTION getNodes1DCart - PURE FUNCTION intersection1DCart(self, r0, v0) RESULT(r) + PURE FUNCTION intersection1DCart(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r r = (/ self%x, 0.D0, 0.D0 /) @@ -383,7 +382,7 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm SUBROUTINE scatterSegm(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index d0f8311..08609cd 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -7,6 +7,9 @@ MODULE moduleMesh1DRad 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 @@ -174,11 +177,11 @@ MODULE moduleMesh1DRad END FUNCTION getNodes1DRad - PURE FUNCTION intersection1DRad(self, r0, v0) RESULT(r) + PURE FUNCTION intersection1DRad(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r r = (/ self%r, 0.D0, 0.D0 /) @@ -317,19 +320,26 @@ MODULE moduleMesh1DRad REAL(8):: ke(1:2,1:2) REAL(8):: Xii(1:3) REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ + REAL(8):: invJ(1), detJ REAL(8):: r, fPsi(1:2) + INTEGER:: l ke = 0.D0 Xii = 0.D0 - dPsi = self%dPsi(Xii) - invJ = self%invJac(Xii, dPsi) - fPsi = self%fPsi(Xii) - r = DOT_PRODUCT(fPsi, self%r) - ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) - ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) - ke = 2.D0*ke*invJ - ke = ke*r*PI2 + 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) + ke = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + r*wSeg(l)/detJ + + END DO + + ke = ke*PI2 END FUNCTION elemKRad @@ -341,17 +351,23 @@ MODULE moduleMesh1DRad REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:2) - REAL(8):: r - REAL(8):: detJ + REAL(8):: detJ, f, r REAL(8):: Xii(1:3) + INTEGER:: l - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) - r = DOT_PRODUCT(fPsi,self%r) ALLOCATE(localF(1:2)) - localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ - localF = localF*r*PI2 + 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 @@ -377,7 +393,7 @@ MODULE moduleMesh1DRad END FUNCTION insideRad SUBROUTINE scatterRad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 40c218c..10b12fe 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -190,9 +190,9 @@ MODULE moduleMesh2DCart self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) !Normal vector - self%normal = (/ self%y(2)-self%y(1), & - self%x(2)-self%x(1), & - 0.D0 /) + self%normal = (/ -(self%y(2)-self%y(1)), & + self%x(2)-self%x(1) , & + 0.D0 /) self%normal = self%normal/NORM2(self%normal) !Boundary index @@ -243,20 +243,21 @@ MODULE moduleMesh2DCart END FUNCTION getNodes2DCart - PURE FUNCTION intersection2DCartEdge(self, r0, v0) RESULT(r) + PURE FUNCTION intersection2DCartEdge(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r - REAL(8), DIMENSION(1:3):: rS !base point of surface - REAL(8):: d + REAL(8), DIMENSION(1:3):: edge0, edgeV + REAL(8):: tI - rS = (/ self%x(1), self%y(1), 0.D0 /) + edge0 = (/self%x(1), self%y(1), 0.D0 /) + edgeV = (/self%x(2), self%y(2), 0.D0 /) - edge0 - d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal) + tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV) - r = r0 + v0*d + r = edge0 + tI*edgeV END FUNCTION intersection2DCartEdge @@ -496,7 +497,7 @@ MODULE moduleMesh2DCart !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -855,7 +856,7 @@ MODULE moduleMesh2DCart !Scatter properties of particles into element SUBROUTINE scatterTria(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 index b0f2ee3..96664df 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 @@ -437,15 +437,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -456,15 +456,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -475,15 +475,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n4%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -494,15 +494,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n4%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e4 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n4%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e4 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -521,16 +521,16 @@ MODULE moduleMesh2DCartRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal + END IF END IF @@ -540,15 +540,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -559,15 +559,15 @@ MODULE moduleMesh2DCartRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index e0f3a27..d48af83 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -191,9 +191,9 @@ MODULE moduleMesh2DCyl self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) !Normal vector - self%normal = (/ self%r(2)-self%r(1), & - self%z(2)-self%z(1), & - 0.D0 /) + self%normal = (/ -(self%r(2)-self%r(1)), & + self%z(2)-self%z(1) , & + 0.D0 /) self%normal = self%normal/NORM2(self%normal) !Boundary index self%boundary => boundary(bt) @@ -221,20 +221,22 @@ MODULE moduleMesh2DCyl END FUNCTION getNodes2DCyl - PURE FUNCTION intersection2DCylEdge(self, r0, v0) RESULT(r) + PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r) + USE moduleMath IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), DIMENSION(1:3), INTENT(in):: r0 REAL(8), DIMENSION(1:3):: r - REAL(8), DIMENSION(1:3):: rS !base point of surface - REAL(8):: d + REAL(8), DIMENSION(1:3):: edge0, edgeV + REAL(8):: tI - rS = (/ self%z(1), self%r(1), 0.D0 /) + edge0 = (/self%z(1), self%r(1), 0.D0 /) + edgeV = (/self%z(2), self%r(2), 0.D0 /) - edge0 - d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal) + tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV) - r = r0 + v0*d + r = edge0 + tI*edgeV END FUNCTION intersection2DCylEdge @@ -507,7 +509,7 @@ MODULE moduleMesh2DCyl !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -874,7 +876,7 @@ MODULE moduleMesh2DCyl !Scatter properties of particles into element SUBROUTINE scatterTria(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 index 379a0fb..cbe3faf 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 @@ -217,6 +217,29 @@ MODULE moduleMesh2DCylRead END SUBROUTINE connectedVolEdge + PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident) + IMPLICIT NONE + + INTEGER, DIMENSION(1:2), INTENT(in):: nodesA, nodesB + LOGICAL:: coincident + INTEGER:: i + + coincident = .FALSE. + DO i = 1, 2 + IF (ANY(nodesA(i) == nodesB)) THEN + coincident = .TRUE. + + ELSE + coincident = .FALSE. + EXIT + + END IF + + END DO + + END FUNCTION coincidentNodes + + SUBROUTINE connectedQuadQuad(elemA, elemB) IMPLICIT NONE @@ -437,15 +460,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -456,15 +479,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -475,15 +498,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n4%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n4%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -494,12 +517,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n4%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e4 => elemB - elemB%e2 => elemA + elemB%e1 => elemA ELSEIF (elemA%n4%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e4 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -518,15 +544,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n1%n == elemB%n1%n .AND. & elemA%n2%n == elemB%n2%n) THEN elemA%e1 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n1%n == elemB%n2%n .AND. & elemA%n2%n == elemB%n1%n) THEN elemA%e1 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -537,15 +563,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n2%n == elemB%n1%n .AND. & elemA%n3%n == elemB%n2%n) THEN elemA%e2 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n2%n == elemB%n2%n .AND. & elemA%n3%n == elemB%n1%n) THEN elemA%e2 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF @@ -556,15 +582,15 @@ MODULE moduleMesh2DCylRead IF (elemA%n3%n == elemB%n1%n .AND. & elemA%n1%n == elemB%n2%n) THEN elemA%e3 => elemB - elemB%e2 => elemA - - !Revers the normal to point inside the domain - elemB%normal = - elemB%normal + elemB%e1 => elemA ELSEIF (elemA%n3%n == elemB%n2%n .AND. & elemA%n1%n == elemB%n1%n) THEN elemA%e3 => elemB - elemB%e1 => elemA + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = - elemB%normal END IF diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 315cb42..ec7f2b1 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -134,6 +134,7 @@ MODULE moduleMesh3DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleMath IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(out):: self @@ -142,6 +143,7 @@ MODULE moduleMesh3DCart INTEGER, INTENT(in):: bt INTEGER, INTENT(in):: physicalSurface REAL(8), DIMENSION(1:3):: r1, r2, r3 + REAL(8), DIMENSION(1:3):: vec1, vec2 INTEGER:: s self%n = n @@ -156,10 +158,14 @@ MODULE moduleMesh3DCart self%y = (/r1(2), r2(2), r3(2)/) self%z = (/r1(3), r2(3), r3(3)/) !Normal vector - self%normal = (/ (self%y(2)-self%y(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)), & - (self%x(2)-self%x(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%x(3)-self%x(1)), & - (self%x(2)-self%x(1))*(self%y(3)-self%y(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)) /) - self%normal = self%normal/NORM2(self%normal) + vec1 = (/ self%x(2) - self%x(1), & + self%y(2) - self%y(1), & + self%z(2) - self%z(1) /) + vec2 = (/ self%x(3) - self%x(1), & + self%y(3) - self%y(1), & + self%z(3) - self%z(1) /) + self%normal = crossProduct(vec1, vec2) + self%normal = normalize(self%normal) !Boundary index self%boundary => boundary(bt) @@ -187,20 +193,21 @@ MODULE moduleMesh3DCart END FUNCTION getNodes3DCartTria - PURE FUNCTION intersection3DCartTria(self, r0, v0) RESULT(r) + PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r) IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0 + REAL(8), INTENT(in):: r0(1:3) REAL(8), DIMENSION(1:3):: r - REAL(8), DIMENSION(1:3):: rS !base point of surface - REAL(8):: d + REAL(8), DIMENSION(1:3):: edge0, edgeV + REAL(8):: tI - rS = (/ self%x(1), self%y(1), self%z(1) /) + edge0 = (/self%x(1), self%y(1), self%z(1) /) + edgeV = (/self%x(2), self%y(2), self%z(2) /) - edge0 - d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal) + tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV) - r = r0 + v0*d + r = edge0 + tI*edgeV END FUNCTION intersection3DCartTria @@ -476,7 +483,7 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra SUBROUTINE scatterTetra(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -577,7 +584,7 @@ MODULE moduleMesh3DCart INTEGER:: nextInt !TODO: Review when connectivity - xiArray = (/ xi(3), xi(2), 1.D0 - xi(1) - xi(2) - xi(3), xi(1) /) + xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /) nextInt = MINLOC(xiArray, 1) NULLIFY(nextElement) SELECT CASE(nextInt) @@ -647,9 +654,11 @@ MODULE moduleMesh3DCart invJ(2,2) = (dx(1)*dz(3) - dx(3)*dz(1)) invJ(2,3) = -(dx(1)*dz(2) - dx(2)*dz(1)) - invJ(3,1) = -(dx(2)*dy(3) - dx(3)*dy(2)) - invJ(3,2) = (dx(1)*dy(3) - dx(3)*dy(1)) - invJ(3,3) = -(dx(1)*dy(2) - dx(2)*dy(1)) + invJ(3,1) = (dx(2)*dy(3) - dx(3)*dy(2)) + invJ(3,2) = -(dx(1)*dy(3) - dx(3)*dy(1)) + invJ(3,3) = (dx(1)*dy(2) - dx(2)*dy(1)) + + invJ = TRANSPOSE(invJ) END FUNCTION invJ3DCart diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 index c6b021a..01aef8e 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 @@ -269,25 +269,25 @@ MODULE moduleMesh3DCartRead !Check surface 2 IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e2 => elemB elemB%e1 => elemA - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e2 => elemA - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e2 => elemB elemB%e3 => elemA - ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e2 => elemB @@ -299,25 +299,25 @@ MODULE moduleMesh3DCartRead !Check surface 3 IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN elemA%e3 => elemB elemB%e1 => elemA - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e2 => elemA - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e3 => elemA - ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB @@ -328,7 +328,7 @@ MODULE moduleMesh3DCartRead END IF !Check surface 4 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (.NOT. ASSOCIATED(elemA%e4)) THEN IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), & (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN @@ -360,77 +360,78 @@ MODULE moduleMesh3DCartRead END SUBROUTINE connectedTetraTetra SUBROUTINE connectedTetraEdge(elemA, elemB) + USE moduleMath IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB + INTEGER:: nodesEdge(1:3) + REAL(8), DIMENSION(1:3):: vec1, vec2 + REAL(8):: normVol(1:3) + + nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /) !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 + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), & + nodesEdge)) THEN - elemA%e1 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(2) - elemA%x(1), & + elemA%y(2) - elemA%y(1), & + elemA%z(2) - elemA%z(1) /) + vec2 = (/ elemA%x(3) - elemA%x(1), & + elemA%y(3) - elemA%y(1), & + elemA%z(3) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - 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 + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e1 => elemB - elemB%e2 => elemA + elemA%e1 => elemB + elemB%e1 => elemA - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal + ELSE + + elemA%e1 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF 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 + IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + nodesEdge)) THEN - elemA%e2 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(3) - elemA%x(2), & + elemA%y(3) - elemA%y(2), & + elemA%z(3) - elemA%z(2) /) + vec2 = (/ elemA%x(4) - elemA%x(2), & + elemA%y(4) - elemA%y(2), & + elemA%z(4) - elemA%z(2) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - 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 + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e2 => elemB - elemB%e2 => elemA + elemA%e2 => elemB + elemB%e1 => elemA - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal + ELSE + + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF END IF @@ -438,66 +439,66 @@ MODULE moduleMesh3DCartRead !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 + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), & + nodesEdge)) THEN - elemA%e3 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(2) - elemA%x(1), & + elemA%y(2) - elemA%y(1), & + elemA%z(2) - elemA%z(1) /) + vec2 = (/ elemA%x(4) - elemA%x(1), & + elemA%y(4) - elemA%y(1), & + elemA%z(4) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - 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 + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e3 => elemB - elemB%e2 => elemA + elemA%e3 => elemB + elemB%e1 => elemA - !Revers the normal to point inside the domain - elemB%normal = -elemB%normal + ELSE + + elemA%e3 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF 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 + IF (.NOT. ASSOCIATED(elemA%e4)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), & + nodesEdge)) THEN - elemA%e4 => elemB - elemB%e1 => elemA + vec1 = (/ elemA%x(3) - elemA%x(1), & + elemA%y(3) - elemA%y(1), & + elemA%z(3) - elemA%z(1) /) + vec2 = (/ elemA%x(4) - elemA%x(1), & + elemA%y(4) - elemA%y(1), & + elemA%z(4) - elemA%z(1) /) + normVol = crossProduct(vec1, vec2) + normVol = normalize(normVol) - 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 + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN - elemA%e4 => elemB - elemB%e2 => elemA + elemA%e4 => elemB + elemB%e1 => elemA + + + ELSE + + elemA%e4 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index f5f7b47..137b256 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -97,10 +97,10 @@ MODULE moduleMesh 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) + PURE FUNCTION intersectionEdge_interface(self, r0) RESULT(r) IMPORT:: meshEdge CLASS(meshEdge), INTENT(in):: self - REAL(8), INTENT(in), DIMENSION(1:3):: r0, v0 + REAL(8), INTENT(in), DIMENSION(1:3):: r0 REAL(8):: r(1:3) END FUNCTION intersectionEdge_interface diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 2252b84..7f5fb88 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -13,30 +13,26 @@ MODULE moduleMeshBoundary !rp = intersection between particle and edge !rpp = final position of particle !vpp = final velocity of particle - REAL(8), DIMENSION(1:3):: rp, rpp, vpp + REAL(8), DIMENSION(1:3):: rp, vpp + REAL(8):: tI REAL(8):: taup !time step for reflecting process !Reflect particle velocity vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal + part%v = vpp - !Computes the intersection between particle and surface - rp = edge%intersection(part%r, part%v) + rp = edge%intersection(part%r) - !Computes the reflection time step - taup = NORM2(part%r - rp)*tau(part%sp) + part%r = 2.D0*(rp - part%r) + part%r - !New position of particle - rpp = rp + vpp*taup - - !assign new parameters to particle - part%r = rpp - part%v = vpp + !particle is assumed to be inside part%n_in = .TRUE. END SUBROUTINE reflection !Absoption in a surface SUBROUTINE absorption(edge, part) + USE moduleCaseParam USE moduleSpecies IMPLICIT NONE @@ -45,7 +41,7 @@ MODULE moduleMeshBoundary REAL(8):: rpp(1:3) !Position of particle projected to the edge REAL(8):: d !Distance from particle to edge - rpp = edge%intersection(part%r, part%v) + rpp = edge%intersection(part%r) d = NORM2(rpp - part%r) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index df9df48..9d50abc 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -200,6 +200,7 @@ MODULE moduleInject USE moduleSolver USE moduleMesh USE moduleRandom + USE moduleErrors IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self @@ -241,7 +242,7 @@ MODULE moduleInject partInj(n)%vol = randomEdge%e2%n ELSE - PRINT *, "ERROR NO VOL ASSOCIATED TO EDGE" + CALL criticalError("No Volume associated to edge", 'addParticles') END IF diff --git a/src/modules/moduleMath.f90 b/src/modules/moduleMath.f90 new file mode 100644 index 0000000..213afa5 --- /dev/null +++ b/src/modules/moduleMath.f90 @@ -0,0 +1,42 @@ +MODULE moduleMath + IMPLICIT NONE + + CONTAINS + !Outer product of two tensors + PURE FUNCTION outerProduct(a,b) RESULT(s) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3), INTENT(in):: a, b + REAL(8):: s(1:3,1:3) + + s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3) + + END FUNCTION outerProduct + + !Cross product of two 3D vectors + PURE FUNCTION crossProduct(a, b) RESULT(c) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3), INTENT(in):: a, b + REAL(8), DIMENSION(1:3):: c + + c = 0.D0 + + c(1) = a(2)*b(3) - a(3)*b(2) + c(2) = -(a(1)*b(3) - a(3)*b(1)) + c(3) = a(1)*b(2) - a(2)*b(1) + + END FUNCTION crossProduct + + !Normalizes a 3D vector + PURE FUNCTION normalize(a) RESULT(b) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3), INTENT(in):: a + REAL(8), DIMENSION(1:3):: b + + b = a / NORM2(a) + + END FUNCTION normalize + +END MODULE moduleMath diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90 index b4b8317..963302d 100644 --- a/src/modules/moduleOutput.f90 +++ b/src/modules/moduleOutput.f90 @@ -31,16 +31,6 @@ MODULE moduleOutput LOGICAL:: emOutput = .FALSE. CONTAINS - FUNCTION outerProduct(a,b) RESULT(s) - IMPLICIT NONE - - REAL(8), DIMENSION(1:3):: a, b - REAL(8):: s(1:3,1:3) - - s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3) - - END FUNCTION outerProduct - FUNCTION tensorTrace(a) RESULT(t) IMPLICIT NONE @@ -56,6 +46,7 @@ MODULE moduleOutput USE moduleConstParam USE moduleRefParam USE moduleSpecies + USE moduleMath IMPLICIT NONE TYPE(outputNode), INTENT(in):: rawValues