diff --git a/src/makefile b/src/makefile index e3801fd..b398134 100644 --- a/src/makefile +++ b/src/makefile @@ -3,7 +3,8 @@ 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 \ $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.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 d45c7a0..020b971 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 @@ -145,30 +148,13 @@ 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) 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 @@ -189,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 /) @@ -322,22 +308,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 +341,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 @@ -381,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/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/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 1d2bdba..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 @@ -153,23 +156,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 @@ -190,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 /) @@ -333,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 @@ -357,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 @@ -393,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/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/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index ac22b28..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 @@ -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 @@ -259,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 @@ -512,7 +497,7 @@ MODULE moduleMesh2DCart !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -871,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 b91627a..96664df 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 @@ -437,12 +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 + 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 @@ -453,12 +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 + 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 @@ -469,12 +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 + 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 @@ -485,12 +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 + 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 @@ -509,13 +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 + 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 @@ -525,12 +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 + 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 @@ -541,12 +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 + 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 8ef5c6d..d48af83 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -191,35 +191,16 @@ 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) 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 @@ -240,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 @@ -459,7 +442,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 @@ -524,7 +509,7 @@ MODULE moduleMesh2DCyl !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -891,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 085c3c3..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,12 +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 + 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 @@ -453,12 +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 + 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 @@ -469,12 +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 + 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 @@ -485,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 @@ -509,12 +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 + 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 @@ -525,12 +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 + 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 @@ -541,12 +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 + 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 1c713c9..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,11 +143,12 @@ 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 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() @@ -156,33 +158,21 @@ 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) 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 @@ -203,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 @@ -492,7 +483,7 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra SUBROUTINE scatterTetra(self, part) - USE moduleOutput + USE moduleMath USE moduleSpecies IMPLICIT NONE @@ -593,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) @@ -663,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 a2197a4..01aef8e 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 @@ -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%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 ((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%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 ((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%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 ((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%n2%n, elemA%n3%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%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 ((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%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 ((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%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 ((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%n1%n, elemA%n2%n, elemA%n4%n/), & + (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN elemA%e3 => elemB elemB%e4 => elemA @@ -388,56 +326,29 @@ MODULE moduleMesh3DCartRead 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/), & + (/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 @@ -449,10 +360,149 @@ 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 (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), & + nodesEdge)) THEN + + 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) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e1 => elemB + elemB%e1 => elemA + + 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 (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), & + nodesEdge)) THEN + + 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) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e2 => elemB + elemB%e1 => elemA + + ELSE + + elemA%e2 => elemB + elemB%e2 => elemA + + !Revers the normal to point inside the domain + elemB%normal = -elemB%normal + + END IF + + END IF + + END IF + + !Check surface 3 + IF (.NOT. ASSOCIATED(elemA%e3)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), & + nodesEdge)) THEN + + 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) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + elemA%e3 => elemB + elemB%e1 => elemA + + 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%e4)) THEN + IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), & + nodesEdge)) THEN + + 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) + + IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN + + 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 + + END IF END SUBROUTINE connectedTetraEdge @@ -465,6 +515,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/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 97602f4..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) @@ -108,6 +104,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 +193,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..9d50abc 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') @@ -196,6 +200,7 @@ MODULE moduleInject USE moduleSolver USE moduleMesh USE moduleRandom + USE moduleErrors IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self @@ -230,12 +235,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 + CALL criticalError("No Volume associated to edge", 'addParticles') + END IF partInj(n)%v = (/ self%v(1)%obj%randomVel(), & diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 13d43e3..96debdd 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') @@ -653,6 +679,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 +697,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/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/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 diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 4b0124e..2927e21 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 @@ -162,6 +169,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 @@ -437,7 +506,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 +527,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 +593,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