diff --git a/src/modules/moduleEM.f95 b/src/modules/moduleEM.f95 index db37941..141b1a8 100644 --- a/src/modules/moduleEM.f95 +++ b/src/modules/moduleEM.f95 @@ -56,7 +56,7 @@ MODULE moduleEM elField = 0.D0 - xi = part%xLog + xi = part%xi elField = mesh%vols(part%e_p)%obj%gatherEF(xi) diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95 index e489a75..4f20b1f 100644 --- a/src/modules/moduleInput.f95 +++ b/src/modules/moduleInput.f95 @@ -86,8 +86,8 @@ MODULE moduleInput v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity !TODO: Make this solver dependent IF (found_r) THEN - L_ref = 1.D0/(sigma_ref*n_ref) !mean free path sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section + L_ref = 1.D0/(sigma_ref*n_ref) !mean free path ELSE L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length @@ -483,7 +483,6 @@ MODULE moduleInput CALL inject(i)%init(i, v, normal, T, flow, units, pt, physicalSurface) END DO - PRINT *, inject(:)%nParticles !Allocate array for injected particles IF (nPartInj > 0) THEN diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95 index 0b6e37f..65389a7 100644 --- a/src/modules/moduleMesh.f95 +++ b/src/modules/moduleMesh.f95 @@ -125,14 +125,17 @@ MODULE moduleMesh !Total weight of particles inside cell REAL(8):: totalWeight = 0.D0 CONTAINS - PROCEDURE(initVol_interface), DEFERRED, PASS:: init - PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter - PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF - PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes - PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF - PROCEDURE(collision_interface), DEFERRED, PASS:: collision - PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell - PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput + PROCEDURE(initVol_interface), DEFERRED, PASS:: init + PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter + PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF + PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes + PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF + PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell + PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log + PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside + PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement + PROCEDURE(collision_interface), DEFERRED, PASS:: collision + PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput END TYPE meshVol @@ -164,7 +167,6 @@ MODULE moduleMesh END FUNCTION gatherEF_interface PURE FUNCTION getNodesVol_interface(self) RESULT(n) - IMPORT:: meshVol CLASS(meshVol), INTENT(in):: self INTEGER, ALLOCATABLE:: n(:) @@ -172,7 +174,6 @@ MODULE moduleMesh END FUNCTION getNodesVol_interface PURE FUNCTION elemF_interface(self, source) RESULT(localF) - IMPORT:: meshVol CLASS(meshVol), INTENT(in):: self REAL(8), INTENT(in):: source(1:) @@ -180,11 +181,13 @@ MODULE moduleMesh END FUNCTION elemF_interface - SUBROUTINE collision_interface(self) + SUBROUTINE nextElement_interface(self, xi, nextElement) IMPORT:: meshVol - CLASS(meshVol), INTENT(inout):: self + CLASS(meshVol), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + CLASS(*), POINTER, INTENT(out):: nextElement - END SUBROUTINE collision_interface + END SUBROUTINE nextElement_interface SUBROUTINE findCell_interface(self, part, oldCell) USE moduleSpecies @@ -196,6 +199,27 @@ MODULE moduleMesh END SUBROUTINE findCell_interface + PURE FUNCTION phy2log_interface(self,r) RESULT(xN) + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: xN(1:3) + + END FUNCTION phy2log_interface + + PURE FUNCTION inside_interface(xi) RESULT(ins) + IMPORT:: meshVol + REAL(8), INTENT(in):: xi(1:3) + LOGICAL:: ins + + END FUNCTION inside_interface + + SUBROUTINE collision_interface(self) + IMPORT:: meshVol + CLASS(meshVol), INTENT(inout):: self + + END SUBROUTINE collision_interface + SUBROUTINE resetOutput_interface(self) IMPORT:: meshVol CLASS(meshVol), INTENT(inout):: self diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95 index ac912c8..69895db 100644 --- a/src/modules/moduleMeshCyl.f95 +++ b/src/modules/moduleMeshCyl.f95 @@ -6,9 +6,13 @@ MODULE moduleMeshCyl USE moduleMesh IMPLICIT NONE - !TODO: Move this, this is horrible - REAL(8), PARAMETER:: w(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) - REAL(8), PARAMETER:: cor(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) + !Values for Gauss integral + REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) + REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) + + REAL(8), PARAMETER:: xi1Tria(1:3) = (/ 1.D0/6.D0, 2.D0/3.D0, 1.D0/6.D0 /) + REAL(8), PARAMETER:: xi2Tria(1:3) = (/ 1.D0/6.D0, 1.D0/6.D0, 2.D0/3.D0 /) + REAL(8), PARAMETER:: wTria(1:3) = (/ 1.D0/6.D0, 1.D0/6.D0, 1.D0/6.D0 /) TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl !Element coordinates @@ -34,9 +38,38 @@ MODULE moduleMeshCyl TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl CONTAINS PROCEDURE, PASS:: collision => collision2DCyl + PROCEDURE, PASS:: findCell => findCellCyl + PROCEDURE, PASS:: detJac => detJCyl + PROCEDURE, PASS:: invJac => invJCyl + PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi + PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi + PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer END TYPE meshVolCyl + ABSTRACT INTERFACE + PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + END FUNCTION fPsi_interface + + PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: dPsi(:,:) + + END FUNCTION dPsi_interface + + PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr) + IMPORT meshVolCyl + CLASS(meshVolCyl), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + + END SUBROUTINE partialDer_interface + + END INTERFACE + TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad !Element coordinates REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 @@ -48,24 +81,58 @@ MODULE moduleMeshCyl CONTAINS PROCEDURE, PASS:: init => initVolQuadCyl + PROCEDURE, PASS:: area => areaQuad + PROCEDURE, NOPASS:: fPsi => fPsiQuad + PROCEDURE, NOPASS:: dPsi => dPsiQuad + PROCEDURE, NOPASS:: dPsiXi1 => dPsiQuadXi1 + PROCEDURE, NOPASS:: dPsiXi2 => dPsiQuadXi2 + PROCEDURE, PASS:: partialDer => partialDerQuad PROCEDURE, PASS:: elemK => elemKQuad PROCEDURE, PASS:: elemF => elemFQuad - PROCEDURE, PASS:: detJac => detJQuad - PROCEDURE, PASS:: invJ => invJQuad - PROCEDURE, PASS:: area => areaQuad PROCEDURE, NOPASS:: weight => weightQuad PROCEDURE, NOPASS:: inside => insideQuad PROCEDURE, PASS:: scatter => scatterQuad PROCEDURE, PASS:: gatherEF => gatherEFQuad PROCEDURE, PASS:: getNodes => getNodesQuad PROCEDURE, PASS:: phy2log => phy2logQuad - PROCEDURE, PASS:: findCell => findCellCylQuad + PROCEDURE, PASS:: nextElement => nextElementQuad PROCEDURE, PASS:: resetOutput => resetOutputQuad END TYPE meshVolCylQuad + TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylTria + !Element coordinates + REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() + !Connectivity to adjacent elements + CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() + REAL(8):: arNodes(1:3) = 0.D0 + + CONTAINS + PROCEDURE, PASS:: init => initVolTriaCyl + PROCEDURE, PASS:: area => areaTria + PROCEDURE, NOPASS:: fPsi => fPsiTria + PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 + PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 + PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, PASS:: elemK => elemKTria + PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, NOPASS:: weight => weightTria + PROCEDURE, NOPASS:: inside => insideTria + PROCEDURE, PASS:: scatter => scatterTria + PROCEDURE, PASS:: gatherEF => gatherEFTria + PROCEDURE, PASS:: getNodes => getNodesTria + PROCEDURE, PASS:: phy2log => phy2logTria + PROCEDURE, PASS:: nextElement => nextElementTria + PROCEDURE, PASS:: resetOutput => resetOutputTria + + END TYPE meshVolCylTria + CONTAINS + !NODE FUNCTIONS !Inits node element SUBROUTINE initNodeCyl(self, n, r) USE moduleSpecies @@ -97,7 +164,8 @@ MODULE moduleMeshCyl END FUNCTION getCoordCyl - !Edge functions + !EDGE FUNCTIONS + !Inits edge element SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface) IMPLICIT NONE @@ -127,6 +195,7 @@ MODULE moduleMeshCyl END SUBROUTINE initEdgeCyl + !Get nodes from edge PURE FUNCTION getNodesCyl(self) RESULT(n) IMPLICIT NONE @@ -138,6 +207,7 @@ MODULE moduleMeshCyl END FUNCTION getNodesCyl + !Calculates a random position in edge FUNCTION randPosCyl(self) RESULT(r) CLASS(meshEdgeCyl), INTENT(in):: self REAL(8):: rnd @@ -153,8 +223,9 @@ MODULE moduleMeshCyl END FUNCTION randPosCyl - !Vol functions - !Quad Volume + !VOLUME FUNCTIONS + !QUAD FUNCTIONS + !Inits quadrilateral element SUBROUTINE initVolQuadCyl(self, n, p) USE moduleRefParam IMPLICIT NONE @@ -162,7 +233,7 @@ MODULE moduleMeshCyl CLASS(meshVolCylQuad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) - REAL(8):: r1(1:3), r2(1:3), r3(1:3), r4(1:3) + REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 self%n = n self%n1 => mesh%nodes(p(1))%obj @@ -190,79 +261,102 @@ MODULE moduleMeshCyl END SUBROUTINE initVolQuadCyl - !TODO: MAKE ELEMENT DEPENDENT - PURE FUNCTION fpsi(xi1,xi2) RESULT(psi) + !Computes element area + PURE SUBROUTINE areaQuad(self) + USE moduleConstParam IMPLICIT NONE - REAL(8),INTENT(in):: xi1,xi2 - REAL(8):: psi(1:4) - psi(1) = (1.D0-xi1)*(1.D0-xi2) - psi(2) = (1.D0+xi1)*(1.D0-xi2) - psi(3) = (1.D0+xi1)*(1.D0+xi2) - psi(4) = (1.D0-xi1)*(1.D0+xi2) - psi = psi*0.25D0 + CLASS(meshVolCylQuad), INTENT(inout):: self + REAL(8):: r, xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:4) - END FUNCTION fpsi + self%volume = 0.D0 + self%arNodes = 0.D0 + !2D 1 point Gauss Quad Integral + xi = 0.D0 + detJ = self%detJac(xi)*8.D0*PI !4*2*pi + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + self%volume = r*detJ + self%arNodes = fPsi*r*detJ + + END SUBROUTINE areaQuad + + !Computes element functions in point xi + PURE FUNCTION fPsiQuad(xi) RESULT(fpsi) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + ALLOCATE(fPsi(1:4)) + + fPsi(1) = (1.D0-xi(1))*(1.D0-xi(2)) + fPsi(2) = (1.D0+xi(1))*(1.D0-xi(2)) + fPsi(3) = (1.D0+xi(1))*(1.D0+xi(2)) + fPsi(4) = (1.D0-xi(1))*(1.D0+xi(2)) + fPsi = fPsi*0.25D0 + + END FUNCTION fPsiQuad + + !Derivative element function at coordinates xi + PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: dPsi(:,:) + + ALLOCATE(dPsi(1:2,1:4)) + + dPsi(1,:) = dPsiQuadXi1(xi(2)) + dPsi(2,:) = dPsiQuadXi2(xi(1)) + + END FUNCTION dPsiQuad !Derivative element function (xi1) - !TODO: MAKE ELEMENT DEPENDENT - PURE FUNCTION dpsiXi1(xi2) RESULT(psiXi1) + PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) IMPLICIT NONE REAL(8),INTENT(in):: xi2 - REAL(8):: psiXi1(1:4) - psiXi1(1) = -(1.D0-xi2) - psiXi1(2) = (1.D0-xi2) - psiXi1(3) = (1.D0+xi2) - psiXi1(4) = -(1.D0+xi2) - psiXi1 = psiXi1*0.25D0 + REAL(8):: dPsiXi1(1:4) + + dPsiXi1(1) = -(1.D0-xi2) + dPsiXi1(2) = (1.D0-xi2) + dPsiXi1(3) = (1.D0+xi2) + dPsiXi1(4) = -(1.D0+xi2) + dPsiXi1 = dPsiXi1*0.25D0 - END FUNCTION dpsiXi1 + END FUNCTION dPsiQuadXi1 !Derivative element function (xi2) - !TODO: MAKE ELEMENT DEPENDENT - PURE FUNCTION dpsiXi2(xi1) RESULT(psiXi2) + PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2) IMPLICIT NONE REAL(8),INTENT(in):: xi1 - REAL(8):: psiXi2(1:4) - psiXi2(1) = -(1.D0-xi1) - psiXi2(2) = -(1.D0+xi1) - psiXi2(3) = (1.D0+xi1) - psiXi2(4) = (1.D0-xi1) - psiXi2 = psiXi2*0.25D0 + REAL(8):: dPsiXi2(1:4) - END FUNCTION dpsiXi2 + dPsiXi2(1) = -(1.D0-xi1) + dPsiXi2(2) = -(1.D0+xi1) + dPsiXi2(3) = (1.D0+xi1) + dPsiXi2(4) = (1.D0-xi1) + dPsiXi2 = dPsiXi2*0.25D0 - !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logQuad(self,r) RESULT(xN) + END FUNCTION dPsiQuadXi2 + + PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - REAL(8):: xO(1:3), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2) - REAL(8):: conv + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr - !Iterative newton method to transform coordinates - conv=1.D0 - xO=0.D0 + dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) - DO WHILE(conv>1.D-4) - dPsi(1,:) = dpsiXi1(xO(2)) - dPsi(2,:) = dpsiXi2(xO(1)) - j = self%invJ(xO, dPsi) - psi = fpsi(xO(1), xO(2)) - f(1) = DOT_PRODUCT(psi,self%z)-r(1) - f(2) = DOT_PRODUCT(psi,self%r)-r(2) - detJ = self%detJac(xO,dPsi) - xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ - conv=MAXVAL(DABS(xN-xO),1) - xO=xN - - END DO - - END FUNCTION phy2logQuad + END SUBROUTINE partialDerQuad !Computes element local stiffness matrix PURE FUNCTION elemKQuad(self) RESULT(ke) @@ -270,7 +364,8 @@ MODULE moduleMeshCyl IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8):: r, xi(1:3), dPsi(1:2,1:4) + REAL(8):: r, xi(1:3) + REAL(8):: fPsi(1:4), dPsi(1:2,1:4) REAL(8):: ke(1:4,1:4) REAL(8):: j(1:2,1:2), detJ INTEGER:: l, m @@ -278,15 +373,16 @@ MODULE moduleMeshCyl ke=0.D0 !Start 2D Gauss Quad Integral DO l=1, 3 - xi(1) = cor(l) - dPsi(2,:) = dpsiXi2(xi(1)) + xi(2) = corQuad(l) + dPsi(1,:) = self%dPsiXi1(xi(2)) DO m = 1, 3 - xi(2) = cor(m) - dPsi(1,:) = dpsiXi1(xi(2)) + xi(1) = corQuad(m) + dPsi(2,:) = self%dPsiXi2(xi(1)) + fPsi = self%fPsi(xi) detJ = self%detJac(xi,dPsi) - j = self%invJ(xi,dPsi) - r = DOT_PRODUCT(fpsi(xi(1),xi(2)),self%r) - ke = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ + j = self%invJac(xi,dPsi) + r = DOT_PRODUCT(fPsi,self%r) + ke = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*wQuad(l)*wQuad(m)/detJ END DO END DO @@ -302,22 +398,23 @@ MODULE moduleMeshCyl CLASS(meshVolCylQuad), INTENT(in):: self REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: r, xi(1:3), psi(1:4) + REAL(8):: r, xi(1:3) + REAL(8):: fPsi(1:4) REAL(8):: detJ, f INTEGER:: l, m ALLOCATE(localF(1:4)) localF = 0.D0 - xi = 0.D0 + xi = 0.D0 DO l=1, 3 + xi(1) = corQuad(l) DO m = 1, 3 - xi(1) = cor(l) - xi(2) = cor(m) - detJ = self%detJac(xi) - psi = fpsi(xi(1),xi(2)) - r = DOT_PRODUCT(psi,self%r) - f = DOT_PRODUCT(psi,source) - localF = localF + r*f*psi*w(l)*w(m)*detJ + xi(2) = corQuad(m) + detJ = self%detJac(xi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + f = DOT_PRODUCT(fPsi,source) + localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ END DO END DO @@ -325,103 +422,30 @@ MODULE moduleMeshCyl END FUNCTION elemFQuad - !Computes element Jacobian determinant - PURE FUNCTION detJQuad(self,xi,dPsi_in) RESULT(dJ) + !Computes weights in the element nodes + PURE FUNCTION weightQuad(xi) RESULT(w) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:2,1:4) - REAL(8):: dJ - REAL(8):: dPsi(1:2,1:4) - REAL(8):: dz(1:2), dr(1:2) - - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in - - ELSE - dPsi(1,:) = dpsiXi1(xi(2)) - dPsi(2,:) = dpsiXi2(xi(1)) - - END IF - - dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) - dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) - dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) - dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) - dJ = dz(1)*dr(2)-dz(2)*dr(1) - - END FUNCTION detJQuad - - PURE FUNCTION invJQuad(self,xi,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dpsi_in(1:2,1:4) - REAL(8):: dPsi(1:2,1:4) - REAL(8):: dz(1:2), dr(1:2) - REAL(8):: j(1:2,1:2) - - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in - - ELSE - dPsi(1,:) = dpsiXi1(xi(2)) - dPsi(2,:) = dpsiXi2(xi(1)) - - END IF - - dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) - dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) - dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) - dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) - j(1,:) = (/ dr(2), -dz(2) /) - j(2,:) = (/ -dr(1), dz(1) /) - - END FUNCTION invJQuad - - SUBROUTINE areaQuad(self)!Computes element area - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylQuad):: self - REAL(8):: r, xi(1:3) - REAL(8):: detJ, fpsi0(1:4) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !2D 1 point Gauss Quad Integral - xi = 0.D0 - detJ = self%detJac(xi)*8.D0*PI !4*2*pi - fpsi0 = fpsi(xi(1),xi(2)) - r = DOT_PRODUCT(fpsi0,self%r) - self%volume = r*detJ - self%arNodes = fpsi0*r*detJ - - END SUBROUTINE areaQuad - - FUNCTION weightQuad(xi1_p, xi2_p) RESULT(w) !Computes weights in the four vertices. '_p' references particle position in the logical space - IMPLICIT NONE - - REAL(8):: xi1_p, xi2_p REAL(8):: w(1:4) - w = fpsi(xi1_p, xi2_p) + w = fPsiQuad(xi) END FUNCTION weightQuad - FUNCTION insideQuad(xi1_p, xi2_p) RESULT(ins) !Checks if a particle is inside a quad element + !Checks if a particle is inside a quad element + PURE FUNCTION insideQuad(xi) RESULT(ins) IMPLICIT NONE - REAL(8):: xi1_p, xi2_p + REAL(8), INTENT(in):: xi(1:3) LOGICAL:: ins - ins = (xi1_p >= -1.D0 .AND. xi1_p <= 1.D0) .AND. & - (xi2_p >= -1.D0 .AND. xi2_p <= 1.D0) + ins = (xi(1) >= -1.D0 .AND. xi(1) <= 1.D0) .AND. & + (xi(2) >= -1.D0 .AND. xi(2) <= 1.D0) END FUNCTION insideQuad + !Scatter properties of particle into element nodes SUBROUTINE scatterQuad(self, part) USE moduleOutput USE moduleSpecies @@ -433,7 +457,7 @@ MODULE moduleMeshCyl REAL(8):: w_p(1:4) REAL(8):: tensorS(1:3,1:3) - w_p = self%weight(part%xLog(1), part%xLog(2)) + w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) vertex => self%n1%output(part%pt) @@ -458,7 +482,8 @@ MODULE moduleMeshCyl END SUBROUTINE scatterQuad - PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)!Computes electric field inside element + !Gathers the electric field at position xi + PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(in):: self @@ -473,16 +498,16 @@ MODULE moduleMeshCyl self%n3%emData%phi, & self%n4%emData%phi /) - dPsi(1,:) = dpsiXi1(xi(2)) - dPsi(2,:) = dpsiXi2(xi(1)) + dPsi = self%dPsi(xi) detJ = self%detJac(xi,dPsi) - j = self%invJ(xi,dPsi) + j = self%invJac(xi,dPsi) EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ EF(3) = 0.D0 END FUNCTION gatherEFQuad + !Gets nodes from quadrilateral element PURE FUNCTION getNodesQuad(self) RESULT(n) IMPLICIT NONE @@ -494,81 +519,64 @@ MODULE moduleMeshCyl END FUNCTION getNodesQuad - RECURSIVE SUBROUTINE findCellCylQuad(self, part, oldCell) - USE moduleSpecies - USE OMP_LIB + !Transforms physical coordinates to element coordinates + PURE FUNCTION phy2logQuad(self,r) RESULT(xN) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(inout):: self - CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell - CLASS(particle), INTENT(inout), TARGET:: part - REAL(8):: xLog(1:3) - REAL(8):: xLogArray(1:4) - CLASS(*), POINTER:: nextElement + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: xN(1:3) + REAL(8):: xO(1:3), detJ, j(1:2,1:2), f(1:2) + REAL(8):: dPsi(1:2,1:4), fPsi(1:4) + REAL(8):: conv + + !Iterative newton method to transform coordinates + conv=1.D0 + xO=0.D0 + + DO WHILE(conv>1.D-4) + dPsi = self%dPsi(xO) + j = self%invJac(xO, dPsi) + fPsi = self%fPsi(xO) + f(1) = DOT_PRODUCT(fPsi,self%z)-r(1) + f(2) = DOT_PRODUCT(fPsi,self%r)-r(2) + detJ = self%detJac(xO,dPsi) + xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ + conv=MAXVAL(DABS(xN-xO),1) + xO=xN + + END DO + + END FUNCTION phy2logQuad + + !Gets the next element for a logical position xi + SUBROUTINE nextElementQuad(self, xi, nextElement) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + CLASS(*), POINTER, INTENT(out):: nextElement + REAL(8):: xiArray(1:4) INTEGER:: nextInt - xLog = self%phy2log(part%r) - !Checks if particle is inside of current cell - IF (self%inside(xLog(1), xLog(2))) THEN - ! IF (PRESENT(oldCell)) THEN - ! !If oldCell, recalculate particle weight, as particle has entered a new cell. - ! part%weight = part%weight*oldCell%volume/self%volume - ! - ! END IF - part%e_p = self%n - part%xLog = xLog - part%n_in = .TRUE. - !Assign particle to listPart_in - CALL OMP_SET_LOCK(self%lock) - CALL self%listPart_in%add(part) - self%totalWeight = self%totalWeight + part%weight - CALL OMP_UNSET_LOCK(self%lock) + xiArray = (/ -xi(2), xi(1), xi(2), -xi(1) /) + nextInt = MAXLOC(xiArray,1) + !Selects the higher value of directions and searches in that direction + NULLIFY(nextElement) + SELECT CASE (nextInt) + CASE (1) + nextElement => self%e1 + CASE (2) + nextElement => self%e2 + CASE (3) + nextElement => self%e3 + CASE (4) + nextElement => self%e4 + END SELECT - ELSE - !If not, searches for a neighbour and repeats the process. - xLogArray = (/ -xLog(2), xLog(1), xLog(2), -xLog(1) /) - nextInt = MAXLOC(xLogArray,1) - !Selects the higher value of directions and searches in that direction - NULLIFY(nextElement) - SELECT CASE (nextInt) - CASE (1) - nextElement => self%e1 - CASE (2) - nextElement => self%e2 - CASE (3) - nextElement => self%e3 - CASE (4) - nextElement => self%e4 - END SELECT - - !Defines the next step - SELECT TYPE(nextElement) - CLASS IS(meshVolCyl) - !Particle moved to new cell, repeat find procedure - CALL nextElement%findCell(part, self) - - CLASS IS (meshEdgeCyl) - !Particle encountered an edge, execute boundary - CALL nextElement%fBoundary(part) - !If particle is still inside the domain, call findCell - IF (part%n_in) THEN - IF(PRESENT(oldCell)) THEN - CALL self%findCell(part, oldCell) - - ELSE - CALL self%findCell(part) - - END IF - END IF - - CLASS DEFAULT - WRITE(*,*) "ERROR, CHECK findCellCylQuad" - - END SELECT - END IF - - END SUBROUTINE findCellCylQuad + END SUBROUTINE nextElementQuad + !Reset the output of nodes in tria element PURE SUBROUTINE resetOutputQuad(self) USE moduleSpecies USE moduleOutput @@ -598,7 +606,473 @@ MODULE moduleMeshCyl END SUBROUTINE resetOutputQuad + !TRIA ELEMENT + !Init tria element + SUBROUTINE initVolTriaCyl(self, n, p) + USE moduleRefParam + IMPLICIT NONE + CLASS(meshVolCylTria), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + REAL(8), DIMENSION(1:3):: r1, r2, r3 + + self%n = n + self%n1 => mesh%nodes(p(1))%obj + self%n2 => mesh%nodes(p(2))%obj + self%n3 => mesh%nodes(p(3))%obj + !Get element coordinates + r1 = self%n1%getCoordinates() + r2 = self%n2%getCoordinates() + r3 = self%n3%getCoordinates() + self%z = (/r1(1), r2(1), r3(1)/) + self%r = (/r1(2), r2(2), r3(2)/) + !Assign node volume + CALL self%area() + self%n1%v = self%n1%v + self%arNodes(1) + self%n2%v = self%n2%v + self%arNodes(2) + self%n3%v = self%n3%v + self%arNodes(3) + + self%sigmaVrelMax = sigma_ref/L_ref**2 + + CALL OMP_INIT_LOCK(self%lock) + + END SUBROUTINE initVolTriaCyl + + !Calculates area for triangular element + PURE SUBROUTINE areaTria(self) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(inout):: self + REAL(8):: r, xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:3) + + self%volume = 0.D0 + self%arNodes = 0.D0 + !2D 1 point Gauss Quad Integral + xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + detJ = self%detJac(xi)*8.D0*PI + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + self%volume = r*detJ + self%arNodes = fPsi*r*detJ + + END SUBROUTINE areaTria + + !Shape functions for triangular element + PURE FUNCTION fPsiTria(xi) RESULT(fPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + ALLOCATE(fPsi(1:3)) + + fPsi(1) = 1.D0 - xi(1) - xi(2) + fPsi(2) = xi(1) + fPsi(3) = xi(2) + + END FUNCTION fPsiTria + + !Derivative element function at coordinates xi + PURE FUNCTION dPsiTria(xi) RESULT(dPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: dPsi(:,:) + + ALLOCATE(dPsi(1:2,1:3)) + + dPsi(1,:) = dPsiTriaXi1(xi(2)) + dPsi(2,:) = dPsiTriaXi2(xi(1)) + + END FUNCTION dPsiTria + + !Derivative element function (xi1) + PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi2 + REAL(8):: dPsiXi1(1:3) + + dPsiXi1(1) = -xi2 + dPsiXi1(2) = 1.D0 + dPsiXi1(3) = 0.D0 + + END FUNCTION dPsiTriaXi1 + + !Derivative element function (xi1) + PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi1 + REAL(8):: dPsiXi2(1:3) + + dPsiXi2(1) = -xi1 + dPsiXi2(2) = 0.D0 + dPsiXi2(3) = 1.D0 + + END FUNCTION dPsiTriaXi2 + + PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + + dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) + + END SUBROUTINE partialDerTria + + !Computes element local stiffness matrix + PURE FUNCTION elemKTria(self) RESULT(ke) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8):: r, xi(1:3) + REAL(8):: fPsi(1:3), dPsi(1:2,1:3) + REAL(8):: ke(1:3,1:3) + REAL(8):: j(1:2,1:2), detJ + INTEGER:: l + + ke=0.D0 + !Start 2D Gauss Quad Integral + DO l=1, 3 + xi(1) = xi1Tria(l) + xi(2) = xi2Tria(l) + dPsi = self%dPsi(xi) + detJ = self%detJac(xi,dPsi) + j = self%invJac(xi,dPsi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + ke = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*wTria(l)/detJ + + END DO + ke = ke*2.D0*PI + + END FUNCTION elemKTria + + !Computes element local source vector + PURE FUNCTION elemFTria(self, source) RESULT(localF) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + REAL(8):: fPsi(1:3) + REAL(8):: r, xi(1:3) + REAL(8):: detJ, f + INTEGER:: l + + ALLOCATE(localF(1:3)) + localF = 0.D0 + xi = 0.D0 + !Start 2D Gauss Quad Integral + DO l=1, 3 + xi(1) = xi1Tria(l) + xi(2) = xi2Tria(l) + detJ = self%detJac(xi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + f = DOT_PRODUCT(fPsi,source) + localF = localF + r*f*fPsi*wTria(l)*detJ + + END DO + localF = localF*2.D0*PI + + END FUNCTION elemFTria + + !Computes weights in the element nodes + PURE FUNCTION weightTria(xi) RESULT(w) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: w(:) + + w = fPsiTria(xi) + + END FUNCTION weightTria + + PURE FUNCTION insideTria(xi) RESULT(ins) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + LOGICAL:: ins + + ins = xi(1) >= 0.D0 .AND. & + xi(2) >= 0.D0 .AND. & + 1.D0 - xi(1) - xi(2) >= 0.D0 + + END FUNCTION insideTria + + !Scatter properties of particles into element + SUBROUTINE scatterTria(self, part) + USE moduleOutput + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + CLASS(particle), INTENT(in):: part + TYPE(outputNode), POINTER:: vertex + REAL(8):: w_p(1:3) + REAL(8):: tensorS(1:3,1:3) + + w_p = self%weight(part%xi) + tensorS = outerProduct(part%v, part%v) + + vertex => self%n1%output(part%pt) + vertex%den = vertex%den + part%weight*w_p(1) + vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS + + vertex => self%n2%output(part%pt) + vertex%den = vertex%den + part%weight*w_p(2) + vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS + + vertex => self%n3%output(part%pt) + vertex%den = vertex%den + part%weight*w_p(3) + vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS + + END SUBROUTINE scatterTria + + !Gathers the electric field at position xi + PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: dPsi(1:2,1:3) + REAL(8):: j(1:2,1:2), detJ + REAL(8):: phi(1:3) + REAL(8):: EF(1:3) + + phi = (/self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi/) + + dPsi = self%dPsi(xi) + detJ = self%detJac(xi,dPsi) + j = self%invJac(xi,dPsi) + EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ + EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ + EF(3) = 0.D0 + + END FUNCTION gatherEFTria + + !Gets nodes from triangular element + PURE FUNCTION getNodesTria(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:3)) + n = (/self%n1%n, self%n2%n, self%n3%n /) + + END FUNCTION getNodesTria + + !Transforms physical coordinates to element coordinates + PURE FUNCTION phy2logTria(self,r) RESULT(xN) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: xN(1:3) + REAL(8):: xO(1:3), detJ, j(1:2,1:2), f(1:2) + REAL(8):: dPsi(1:2,1:3), fPsi(1:3) + REAL(8):: conv + + !Iterative newton method to transform coordinates + !TODO: Try to find a direct alternative + conv=1.D0 + xO=(/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + + DO WHILE(conv>1.D-4) + dPsi = self%dPsi(xO) + j = self%invJac(xO, dPsi) + fPsi = self%fPsi(xO) + f(1) = DOT_PRODUCT(fPsi,self%z)-r(1) + f(2) = DOT_PRODUCT(fPsi,self%r)-r(2) + detJ = self%detJac(xO,dPsi) + xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ + conv=MAXVAL(DABS(xN-xO),1) + xO=xN + + END DO + + END FUNCTION phy2logTria + + SUBROUTINE nextElementTria(self, xi, nextElement) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + CLASS(*), POINTER, INTENT(out):: nextElement + REAL(8):: xiArray(1:3) + INTEGER:: nextInt + + xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) + nextInt = MINLOC(xiArray,1) + NULLIFY(nextElement) + SELECT CASE (nextInt) + CASE (1) + nextElement => self%e1 + CASE (2) + nextElement => self%e2 + CASE (3) + nextElement => self%e3 + END SELECT + + END SUBROUTINE nextElementTria + + !Reset the output of nodes in tria element + PURE SUBROUTINE resetOutputTria(self) + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(inout):: self + INTEGER:: k + + DO k = 1, nSpecies + self%n1%output(k)%den = 0.D0 + self%n1%output(k)%mom = 0.D0 + self%n1%output(k)%tensorS = 0.D0 + + self%n2%output(k)%den = 0.D0 + self%n2%output(k)%mom = 0.D0 + self%n2%output(k)%tensorS = 0.D0 + + self%n3%output(k)%den = 0.D0 + self%n3%output(k)%mom = 0.D0 + self%n3%output(k)%tensorS = 0.D0 + + END DO + + END SUBROUTINE resetOutputTria + + + !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS + !Computes element Jacobian determinant + PURE FUNCTION detJCyl(self, xi, dPsi_in) RESULT(dJ) + IMPLICIT NONE + + CLASS(meshVolCyl), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(:,:) + REAL(8), ALLOCATABLE:: dPsi(:,:) + REAL(8):: dJ + REAL(8):: dz(1:2), dr(1:2) + + IF(PRESENT(dPsi_in)) THEN + dPsi = dPsi_in + + ELSE + dPsi = self%dPsi(xi) + + END IF + + CALL self%partialDer(dPsi, dz, dr) + dJ = dz(1)*dr(2)-dz(2)*dr(1) + + END FUNCTION detJCyl + + !Computes element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJCyl(self,xi,dPsi_in) RESULT(invJ) + IMPLICIT NONE + + CLASS(meshVolCyl), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(:,:) + REAL(8):: dPsi(1:2,1:4) + REAL(8):: dz(1:2), dr(1:2) + REAL(8):: invJ(1:2,1:2) + + IF(PRESENT(dPsi_in)) THEN + dPsi=dPsi_in + + ELSE + dPsi = self%dPsi(xi) + + END IF + + CALL self%partialDer(dPsi, dz, dr) + invJ(1,:) = (/ dr(2), -dz(2) /) + invJ(2,:) = (/ -dr(1), dz(1) /) + + END FUNCTION invJCyl + + !Find next cell for particle + RECURSIVE SUBROUTINE findCellCyl(self, part, oldCell) + USE moduleSpecies + USE OMP_LIB + IMPLICIT NONE + + CLASS(meshVolCyl), INTENT(inout):: self + CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell + CLASS(particle), INTENT(inout), TARGET:: part + REAL(8):: xi(1:3) + CLASS(*), POINTER:: nextElement + + xi = self%phy2log(part%r) + !Checks if particle is inside of current cell + IF (self%inside(xi)) THEN + ! IF (PRESENT(oldCell)) THEN + ! !If oldCell, recalculate particle weight, as particle has entered a new cell. + ! part%weight = part%weight*oldCell%volume/self%volume + ! + ! END IF + part%e_p = self%n + part%xi = xi + part%n_in = .TRUE. + !Assign particle to listPart_in + CALL OMP_SET_LOCK(self%lock) + CALL self%listPart_in%add(part) + self%totalWeight = self%totalWeight + part%weight + CALL OMP_UNSET_LOCK(self%lock) + + ELSE + !If not, searches for a neighbour and repeats the process. + CALL self%nextElement(xi, nextElement) + !Defines the next step + SELECT TYPE(nextElement) + CLASS IS(meshVolCyl) + !Particle moved to new cell, repeat find procedure + CALL nextElement%findCell(part, self) + + CLASS IS (meshEdgeCyl) + !Particle encountered an edge, execute boundary + CALL nextElement%fBoundary(part) + !If particle is still inside the domain, call findCell + IF (part%n_in) THEN + IF(PRESENT(oldCell)) THEN + CALL self%findCell(part, oldCell) + + ELSE + CALL self%findCell(part) + + END IF + END IF + + CLASS DEFAULT + WRITE(*,*) "ERROR, CHECK findCellCylQuad" + + END SELECT + END IF + + END SUBROUTINE findCellCyl + + !Computes collisions in element SUBROUTINE collision2DCyl(self) USE moduleCollisions USE moduleSpecies diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95 index dce36a4..8ba5509 100644 --- a/src/modules/moduleSpecies.f95 +++ b/src/modules/moduleSpecies.f95 @@ -31,7 +31,7 @@ MODULE moduleSpecies REAL(8):: v(1:3) !Velocity INTEGER:: pt !Particle species id INTEGER:: e_p !Index of element in which the particle is located - REAL(8):: xLog(1:3) !Logical coordinates of particle in element e_p. + REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p. LOGICAL:: n_in !Flag that indicates if a particle is in the domain REAL(8):: weight=0.D0 !weight of particle REAL(8):: qm = 0.D0 !charge over mass