From a2631f6b78213ce390ba1bed5d9fb2fa8cf4e285 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 3 Apr 2021 09:20:46 +0200 Subject: [PATCH] Impliementation of a collision mesh which is independent for the mesh used to scatter particles and compute the EM field. --- runs/Argon_Expansion/CX_case.json | 2 +- runs/cylFlow/input.json | 2 +- src/fpakc.f90 | 12 +- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 21 +- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 21 +- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 34 +- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 34 +- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 25 +- .../mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 | 262 +++++++------- .../inout/gmsh2/moduleMeshOutputGmsh2.f90 | 20 +- src/modules/mesh/moduleMesh.f90 | 330 +++++++++++++----- src/modules/moduleCollisions.f90 | 1 + src/modules/moduleCompTime.f90 | 13 +- src/modules/moduleInject.f90 | 23 ++ src/modules/moduleInput.f90 | 173 +++++---- src/modules/moduleList.f90 | 2 +- src/modules/moduleOutput.f90 | 2 +- src/modules/moduleSolver.f90 | 26 +- src/modules/moduleSpecies.f90 | 1 + 19 files changed, 636 insertions(+), 368 deletions(-) diff --git a/runs/Argon_Expansion/CX_case.json b/runs/Argon_Expansion/CX_case.json index 463c8af..d382a5d 100644 --- a/runs/Argon_Expansion/CX_case.json +++ b/runs/Argon_Expansion/CX_case.json @@ -3,7 +3,7 @@ "path": "./runs/Argon_Expansion/", "triggerOutput": 10, "cpuTime": false, - "numColl": false, + "numColl": true, "folder": "CX_case" }, "geometry": { diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json index b61af96..77b274c 100644 --- a/runs/cylFlow/input.json +++ b/runs/cylFlow/input.json @@ -3,7 +3,7 @@ "path": "./runs/cylFlow/", "triggerOutput": 10, "cpuTime": true, - "numColl": false + "numColl": true }, "geometry": { "type": "2DCyl", diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 3e68a77..4028ba9 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -4,6 +4,7 @@ PROGRAM fpakc USE moduleErrors USE moduleInject USE moduleSolver + USE moduleMesh USE moduleCompTime USE moduleCaseParam USE OMP_LIB @@ -67,10 +68,19 @@ PROGRAM fpakc tColl = omp_get_wtime() !$OMP END SINGLE - CALL doCollisions() + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions() !$OMP SINGLE tColl = omp_get_wtime() - tColl + + !Coulomb scattering + tCoul = omp_get_wTime() + !$OMP END SINGLE + + IF (ASSOCIATED(mesh%doCoulomb)) CALL mesh%doCoulomb() + + !$OMP SINGLE + tCoul = omp_get_wTime() - tCoul !Reset particles tReset = omp_get_wtime() diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 12898d5..a6e4248 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -198,18 +198,19 @@ MODULE moduleMesh1DCart !VOLUME FUNCTIONS !SEGMENT FUNCTIONS !Init segment element - SUBROUTINE initVol1DCartSegm(self, n, p) + SUBROUTINE initVol1DCartSegm(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol1DCartSegm), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2 self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -525,7 +526,7 @@ MODULE moduleMesh1DCart SUBROUTINE connectMesh1DCart(self) IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et DO e = 1, self%numVols @@ -538,11 +539,15 @@ MODULE moduleMesh1DCart END DO - !Connect Vol-Edge - DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + SELECT TYPE(self) + TYPE IS(meshParticles) + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) - END DO + END DO + + END SELECT END DO diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index dc47c62..2bd9d5b 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -200,18 +200,19 @@ MODULE moduleMesh1DRad !VOLUME FUNCTIONS !SEGMENT FUNCTIONS !Init segment element - SUBROUTINE initVol1DRadSegm(self, n, p) + SUBROUTINE initVol1DRadSegm(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol1DRadSegm), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2 self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -536,7 +537,7 @@ MODULE moduleMesh1DRad SUBROUTINE connectMesh1DRad(self) IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et DO e = 1, self%numVols @@ -549,11 +550,15 @@ MODULE moduleMesh1DRad END DO - !Connect Vol-Edge - DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + SELECT TYPE(self) + TYPE IS(meshParticles) + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) - END DO + END DO + + END SELECT END DO diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 11e6454..89704d6 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -283,20 +283,21 @@ MODULE moduleMesh2DCart !VOLUME FUNCTIONS !QUAD FUNCTIONS !Inits quadrilateral element - SUBROUTINE initVolQuad2DCart(self, n, p) + SUBROUTINE initVolQuad2DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj - self%n4 => mesh%nodes(p(4))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj + self%n3 => nodes(p(3))%obj + self%n4 => nodes(p(4))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -631,22 +632,23 @@ MODULE moduleMesh2DCart !TRIA ELEMENT !Init tria element - SUBROUTINE initVolTria2DCart(self, n, p) + SUBROUTINE initVolTria2DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3 !Assign node index self%n = n !Assign nodes to element - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj + self%n3 => nodes(p(3))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -1022,7 +1024,7 @@ MODULE moduleMesh2DCart SUBROUTINE connectMesh2DCart(self) IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et DO e = 1, self%numVols @@ -1035,11 +1037,15 @@ MODULE moduleMesh2DCart END DO - !Connect Vol-Edge - DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + SELECT TYPE(self) + TYPE IS(meshParticles) + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) - END DO + END DO + + END SELECT END DO diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 362380d..2784414 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -271,20 +271,21 @@ MODULE moduleMesh2DCyl !VOLUME FUNCTIONS !QUAD FUNCTIONS !Inits quadrilateral element - SUBROUTINE initVolQuad2DCyl(self, n, p) + SUBROUTINE initVolQuad2DCyl(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol2DCylQuad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj - self%n4 => mesh%nodes(p(4))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj + self%n3 => nodes(p(3))%obj + self%n4 => nodes(p(4))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -652,22 +653,23 @@ MODULE moduleMesh2DCyl !TRIA ELEMENT !Init tria element - SUBROUTINE initVolTria2DCyl(self, n, p) + SUBROUTINE initVolTria2DCyl(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol2DCylTria), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3 !Assign node index self%n = n !Assign nodes to element - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj + self%n3 => nodes(p(3))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -1052,7 +1054,7 @@ MODULE moduleMesh2DCyl SUBROUTINE connectMesh2DCyl(self) IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et DO e = 1, self%numVols @@ -1065,11 +1067,15 @@ MODULE moduleMesh2DCyl END DO - !Connect Vol-Edge - DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + SELECT TYPE(self) + TYPE IS(meshParticles) + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) - END DO + END DO + + END SELECT END DO diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 3e3b970..ad5c65d 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -248,21 +248,22 @@ MODULE moduleMesh3DCart !VOLUME FUNCTIONS !TETRA FUNCTIONS !Inits tetrahedron element - SUBROUTINE initVolTetra3DCart(self, n, p) + SUBROUTINE initVolTetra3DCart(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node REAL(8):: volNodes(1:4) !Volume of each node self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj - self%n4 => mesh%nodes(p(4))%obj + self%n1 => nodes(p(1))%obj + self%n2 => nodes(p(2))%obj + self%n3 => nodes(p(3))%obj + self%n4 => nodes(p(4))%obj !Get element coordinates r1 = self%n1%getCoordinates() r2 = self%n2%getCoordinates() @@ -707,7 +708,7 @@ MODULE moduleMesh3DCart SUBROUTINE connectMesh3DCart(self) IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER:: e, et DO e = 1, self%numVols @@ -720,11 +721,15 @@ MODULE moduleMesh3DCart END DO - !Connect Vol-Edge - DO et = 1, self%numEdges - CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) + SELECT TYPE(self) + TYPE IS(meshParticles) + !Connect Vol-Edge + DO et = 1, self%numEdges + CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj) - END DO + END DO + + END SELECT END DO diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 7cfdb01..59d7064 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -7,11 +7,15 @@ MODULE moduleMeshInputGmsh2 USE moduleMeshOutputGmsh2 IMPLICIT NONE - TYPE(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout), TARGET:: self - self%printOutput => printOutputGmsh2 - self%printColl => printCollGmsh2 - self%printEM => printEMGmsh2 + IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollGmsh2 + SELECT TYPE(self) + TYPE IS(meshParticles) + self%printOutput => printOutputGmsh2 + self%printEM => printEMGmsh2 + + END SELECT self%readMesh => readGmsh2 END SUBROUTINE initGmsh2 @@ -26,12 +30,13 @@ MODULE moduleMeshInputGmsh2 USE moduleBoundary IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: x1, x2, x3 !3 generic coordinates + REAL(8):: r(1:3) !3 generic coordinates INTEGER, ALLOCATABLE:: p(:) !Array for nodes INTEGER:: e = 0, n = 0, eTemp = 0, elemType = 0, bt = 0 INTEGER:: totalNumElem + INTEGER:: numEdges INTEGER:: boundaryType !Read mesh @@ -48,39 +53,44 @@ MODULE moduleMeshInputGmsh2 !Allocate required matrices and vectors ALLOCATE(self%nodes(1:self%numNodes)) - ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes)) - ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes)) - self%K = 0.D0 - self%IPIV = 0 + SELECT TYPE(self) + TYPE IS(meshParticles) + ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes)) + ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes)) + self%K = 0.D0 + self%IPIV = 0 + + END SELECT !Read the nodes information DO e = 1, self%numNodes - READ(10, *) n, x1, x2, x3 + READ(10, *) n, r(1), r(2), r(3) SELECT CASE(self%geometry) CASE("3DCart") ALLOCATE(meshNode3Dcart::self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x1, x2, x3 /)) CASE("2DCyl") ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x1, x2, 0.D0 /)) + r(3) = 0.D0 CASE("2DCart") ALLOCATE(meshNode2DCart:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x1, x2, 0.D0 /)) + r(3) = 0.D0 CASE("1DRad") ALLOCATE(meshNode1DRad:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x1, 0.D0, 0.D0 /)) + r(2:3) = 0.D0 CASE("1DCart") ALLOCATE(meshNode1DCart:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/x1, 0.D0, 0.D0 /)) + r(2:3) = 0.D0 END SELECT + CALL self%nodes(n)%obj%init(n, r) END DO + !Skip comments READ(10, *) READ(10, *) @@ -89,118 +99,116 @@ MODULE moduleMeshInputGmsh2 READ(10, *) totalNumElem !conts edges and volume elements - self%numEdges = 0 - DO e = 1, totalNumElem - READ(10, *) eTemp, elemType - SELECT CASE(self%geometry) - CASE("3DCart") - !Element type 2 is triangle in gmsh - IF (elemType == 2) self%numEdges = e + SELECT TYPE(self) + TYPE IS(meshParticles) + self%numEdges = 0 + DO e = 1, totalNumElem + READ(10, *) eTemp, elemType + SELECT CASE(self%geometry) + CASE("3DCart") + !Element type 2 is triangle in gmsh + IF (elemType == 2) self%numEdges = e - CASE("2DCyl","2DCart") - !Element type 1 is segment in Gmsh - IF (elemType == 1) self%numEdges = e + CASE("2DCyl","2DCart") + !Element type 1 is segment in Gmsh + IF (elemType == 1) self%numEdges = e - CASE("1DRad","1DCart") - !Element type 15 is physical point in Gmsh - IF (elemType == 15) self%numEdges = e - - END SELECT - - END DO - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = TotalnumElem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beggining to read elements - DO e=1, totalNumElem - BACKSPACE(10) - END DO - - !Reads edges - DO e=1, self%numEdges - !Reads the edge according to the geometry - SELECT CASE(self%geometry) - CASE("3DCart") - READ(10, *) n, elemType, eTemp, boundaryType - BACKSPACE(10) - - !Associate boundary condition procedure. - bt = getBoundaryID(boundaryType) - - SELECT CASE(elemType) - CASE(2) - !Triangular surface - ALLOCATE(p(1:3)) - - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3) - - ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:3), bt, boundaryType) - - DEALLOCATE(p) + CASE("1DRad","1DCart") + !Element type 15 is physical point in Gmsh + IF (elemType == 15) self%numEdges = e END SELECT - CASE("2DCyl") - ALLOCATE(p(1:2)) + END DO - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) + !Substract the number of edges to the total number of elements + !to obtain the number of volume elements + self%numVols = TotalnumElem - self%numEdges + ALLOCATE(self%edges(1:self%numEdges)) + numEdges = self%numEdges - ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) + !Go back to the beggining to read elements + DO e=1, totalNumElem + BACKSPACE(10) + END DO - CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) + TYPE IS(meshCollisions) + self%numVols = TotalnumElem + numEdges = 0 + END SELECT + + !Allocates arrays + ALLOCATE(self%vols(1:self%numVols)) + + SELECT TYPE(self) + TYPE IS(meshParticles) + !Reads edges + DO e=1, self%numEdges + !Reads the edge according to the geometry + SELECT CASE(self%geometry) + CASE("3DCart") + READ(10, *) n, elemType, eTemp, boundaryType + BACKSPACE(10) + + !Associate boundary condition procedure. + bt = getBoundaryID(boundaryType) + + SELECT CASE(elemType) + CASE(2) + !Triangular surface + ALLOCATE(p(1:3)) + + READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3) + + ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj) + + END SELECT + + CASE("2DCyl") + ALLOCATE(p(1:2)) + + READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) + !Associate boundary condition procedure. + bt = getBoundaryId(boundaryType) + + ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) + + CASE("2DCart") + ALLOCATE(p(1:2)) + + READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) + !Associate boundary condition procedure. + bt = getBoundaryId(boundaryType) + + ALLOCATE(meshEdge2DCart:: self%edges(e)%obj) + + CASE("1DRad") + ALLOCATE(p(1:1)) + + READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) + !Associate boundary condition + bt = getBoundaryId(boundaryType) + + ALLOCATE(meshEdge1DRad:: self%edges(e)%obj) + + CASE("1DCart") + ALLOCATE(p(1:1)) + + READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) + !Associate boundary condition + bt = getBoundaryId(boundaryType) + + ALLOCATE(meshEdge1DCart:: self%edges(e)%obj) + + END SELECT + + CALL self%edges(e)%obj%init(n, p, bt, boundaryType) DEALLOCATE(p) - CASE("2DCart") - ALLOCATE(p(1:2)) + END DO - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge2DCart:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) - - DEALLOCATE(p) - - CASE("1DRad") - ALLOCATE(p(1:1)) - - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge1DRad:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) - - DEALLOCATE(p) - - CASE("1DCart") - ALLOCATE(p(1:1)) - - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - - ALLOCATE(meshEdge1DCart:: self%edges(e)%obj) - - CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) - - DEALLOCATE(p) - - END SELECT - - END DO + END SELECT !Read and initialize volumes DO e = 1, self%numVols @@ -216,12 +224,9 @@ MODULE moduleMeshInputGmsh2 ALLOCATE(p(1:4)) READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4) ALLOCATE(meshVol3DCartTetra:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) END SELECT - DEALLOCATE(p) - CASE("2DCyl") READ(10,*) n, elemType BACKSPACE(10) @@ -232,19 +237,15 @@ MODULE moduleMeshInputGmsh2 ALLOCATE(p(1:3)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) CASE (3) !Quadrilateral element ALLOCATE(p(1:4)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) END SELECT - DEALLOCATE(p) - CASE("2DCart") READ(10,*) n, elemType BACKSPACE(10) @@ -255,41 +256,36 @@ MODULE moduleMeshInputGmsh2 ALLOCATE(p(1:3)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) CASE (3) !Quadrilateral element ALLOCATE(p(1:4)) READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) END SELECT - DEALLOCATE(p) - CASE("1DRad") ALLOCATE(p(1:2)) READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) ALLOCATE(meshVol1DRadSegm:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) - - DEALLOCATE(p) CASE("1DCart") ALLOCATE(p(1:2)) READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) ALLOCATE(meshVol1DCartSegm:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) - - DEALLOCATE(p) END SELECT + CALL self%vols(e)%obj%init(n - numEdges, p, self%nodes) + DEALLOCATE(p) + END DO + CLOSE(10) + END SUBROUTINE readGmsh2 END MODULE moduleMeshInputGmsh2 diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 902206d..bec5449 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -9,7 +9,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleOutput IMPLICIT NONE - CLASS(meshParticle), INTENT(in):: self + CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t INTEGER:: n, i TYPE(outputFormat):: output(1:self%numNodes) @@ -95,13 +95,25 @@ MODULE moduleMeshOutputGmsh2 USE moduleOutput IMPLICIT NONE - CLASS(meshParticle), INTENT(in):: self + CLASS(meshGeneric), INTENT(in):: self + INTEGER:: numEdges INTEGER, INTENT(in):: t INTEGER:: n REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName CHARACTER (LEN=iterationDigits):: tstring + SELECT TYPE(self) + TYPE IS(meshParticles) + numEdges = self%numEdges + + TYPE IS(meshCollisions) + numEdges = 0 + + CLASS DEFAULT + numEdges = 0 + + END SELECT IF (collOutput) THEN time = DBLE(t)*tauMin*ti_ref @@ -123,7 +135,7 @@ MODULE moduleMeshOutputGmsh2 WRITE(60, *) 1 WRITE(60, *) self%numVols DO n=1, self%numVols - WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl + WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%nColl END DO WRITE(60, "(A)") '$EndElementData' @@ -141,7 +153,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleOutput IMPLICIT NONE - CLASS(meshParticle), INTENT(in):: self + CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t INTEGER:: n, e REAL(8):: time diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index e80b957..93c53d4 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -65,6 +65,8 @@ MODULE moduleMesh TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge !Connectivity to vols CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL() + !Connectivity to vols in meshColl + CLASS(meshVol), POINTER:: eColl => NULL() !Normal vector REAL(8):: normal(1:3) !Weight for random injection of particles @@ -153,8 +155,6 @@ MODULE moduleMesh INTEGER(KIND=OMP_LOCK_KIND):: lock !Number of collisions per volume INTEGER:: nColl = 0 - !Collisional fraction - REAL(8):: collFrac = 0.D0 !Total weight of particles inside cell REAL(8):: totalWeight = 0.D0 CONTAINS @@ -169,16 +169,17 @@ MODULE moduleMesh PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement - PROCEDURE, PASS:: collision END TYPE meshVol ABSTRACT INTERFACE - SUBROUTINE initVol_interface(self, n, p) + SUBROUTINE initVol_interface(self, n, p, nodes) IMPORT:: meshVol + IMPORT meshNodeCont CLASS(meshVol), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) + TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) END SUBROUTINE initVol_interface @@ -260,80 +261,143 @@ MODULE moduleMesh END TYPE meshVolCont - !Particle mesh - TYPE, PUBLIC:: meshParticle - INTEGER:: numEdges, numNodes, numVols - !Array of nodes - TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) - !Array of boundary elements - TYPE(meshEdgeCont), ALLOCATABLE:: edges(:) - !Array of volume elements - TYPE(meshVolCont), ALLOCATABLE:: vols(:) + !Generic mesh type + TYPE, ABSTRACT:: meshGeneric !Geometry of the mesh CHARACTER(:), ALLOCATABLE:: geometry + !Number of elements + INTEGER:: numNodes, numVols + !Array of nodes + TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) + !Array of volume elements + TYPE(meshVolCont), ALLOCATABLE:: vols(:) + PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() + PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() + PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() + CONTAINS + PROCEDURE, PASS:: doCollisions + + END TYPE + + ABSTRACT INTERFACE + !Reads the mesh from a file + SUBROUTINE readMesh_interface(self, filename) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + + END SUBROUTINE readMesh_interface + + !Connects volume and edges to the mesh + SUBROUTINE connectMesh_interface(self) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(inout):: self + + END SUBROUTINE connectMesh_interface + + !Prints number of collisions in each volume + SUBROUTINE printColl_interface(self, t) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + + END SUBROUTINE printColl_interface + + END INTERFACE + + !Particle mesh + TYPE, EXTENDS(meshGeneric), PUBLIC:: meshParticles + INTEGER:: numEdges + !Array of boundary elements + TYPE(meshEdgeCont), ALLOCATABLE:: edges(:) !Global stiffness matrix REAL(8), ALLOCATABLE, DIMENSION(:,:):: K !Permutation matrix for P L U factorization INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() - PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() - PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() - PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() + PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() CONTAINS PROCEDURE, PASS:: constructGlobalK - END TYPE meshParticle + END TYPE meshParticles ABSTRACT INTERFACE + !Perform Coulomb Scattering + SUBROUTINE doCoulomb_interface(self) + IMPORT meshParticles + + CLASS(meshParticles), INTENT(inout):: self + + END SUBROUTINE doCoulomb_interface + !Prints Species data SUBROUTINE printOutput_interface(self, t) - IMPORT meshParticle + IMPORT meshParticles - CLASS(meshParticle), INTENT(in):: self + CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printOutput_interface - !Prints number of collisions - SUBROUTINE printColl_interface(self, t) - IMPORT meshParticle - - CLASS(meshParticle), INTENT(in):: self - INTEGER, INTENT(in):: t - - END SUBROUTINE printColl_interface - !Prints EM info SUBROUTINE printEM_interface(self, t) - IMPORT meshParticle + IMPORT meshParticles - CLASS(meshParticle), INTENT(in):: self + CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface - !Reads the mesh from a file - SUBROUTINE readMesh_interface(self, filename) - IMPORT meshParticle - - CLASS(meshParticle), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - - END SUBROUTINE readMesh_interface - - SUBROUTINE connectMesh_interface(self) - IMPORT meshParticle - - CLASS(meshParticle), INTENT(inout):: self - - END SUBROUTINE connectMesh_interface - END INTERFACE - !Particle mesh - TYPE(meshParticle), TARGET:: mesh + TYPE(meshParticles), TARGET:: mesh + + !Collision (MCC) mesh + TYPE, EXTENDS(meshGeneric):: meshCollisions + CONTAINS + + END TYPE meshCollisions + + TYPE(meshCollisions), TARGET:: meshColl + + ABSTRACT INTERFACE + SUBROUTINE readMeshColl_interface(self, filename) + IMPORT meshCollisions + + CLASS(meshCollisions), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + + END SUBROUTINE readMeshColl_interface + + SUBROUTINE connectMeshColl_interface(self) + IMPORT meshParticles + + CLASS(meshParticles), INTENT(inout):: self + + END SUBROUTINE connectMeshColl_interface + + END INTERFACE + + !Pointer to mesh used for MC collisions + CLASS(meshGeneric), POINTER:: meshForMCC => NULL() + + !Procedure to find a volume for a particle in meshColl + PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL() + + ABSTRACT INTERFACE + SUBROUTINE findCellColl_interface(part) + USE moduleSpecies + + TYPE(particle), INTENT(inout):: part + + END SUBROUTINE findCellColl_interface + + END INTERFACE CONTAINS !Reset the output of node @@ -362,8 +426,8 @@ MODULE moduleMesh IMPLICIT NONE CLASS(meshVol), INTENT(inout):: self - CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell CLASS(particle), INTENT(inout), TARGET:: part + CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell REAL(8):: xi(1:3) CLASS(meshElement), POINTER:: nextElement @@ -408,12 +472,96 @@ MODULE moduleMesh CALL criticalError("No connectivity found for element", "findCell") END SELECT + END IF END SUBROUTINE findCell + !If Coll and Particle are the same, simply copy the part%vol into part%volColl + SUBROUTINE findCellSameMesh(part) + USE moduleSpecies + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + + part%volColl = part%vol + + END SUBROUTINE findCellSameMesh + + !TODO: try to combine this with the findCell for a regular mesh + !Find the volume in which particle reside in the mesh for collisions + SUBROUTINE findCellCollMesh(part) + USE moduleSpecies + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + LOGICAL:: found + CLASS(meshVol), POINTER:: vol + REAL(8), DIMENSION(1:3):: xii + CLASS(meshElement), POINTER:: nextElement + + found = .FALSE. + + vol => meshColl%vols(part%volColl)%obj + DO WHILE(.NOT. found) + xii = vol%phy2log(part%r) + IF (vol%inside(xii)) THEN + part%volColl = vol%n + CALL OMP_SET_LOCK(vol%lock) + CALL vol%listPart_in%add(part) + vol%totalWeight = vol%totalWeight + part%weight + CALL OMP_UNSET_LOCK(vol%lock) + found = .TRUE. + + ELSE + CALL vol%nextElement(xii, nextElement) + SELECT TYPE(nextElement) + CLASS IS(meshVol) + !Try next element + vol => nextElement + + CLASS DEFAULT + !Should never happend, but just in case, stops loops + found = .TRUE. + + END SELECT + + END IF + + END DO + + END SUBROUTINE findCellCollMesh + + !Returns index of volume associated to a position (if any) + !If no voulme is found, returns 0 + !WARNING: This function is slow and should only be used in initialization phase + FUNCTION findCellBrute(self, r) RESULT(nVol) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(in):: self + REAL(8), DIMENSION(1:3), INTENT(in):: r + INTEGER:: nVol + INTEGER:: e + REAL(8), DIMENSION(1:3):: xii + + !Inits RESULT + nVol = 0 + + DO e = 1, self%numVols + xii = self%vols(e)%obj%phy2log(r) + IF(self%vols(e)%obj%inside(xii)) THEN + nVol = self%vols(e)%obj%n + EXIT + + END IF + + END DO + + END FUNCTION findCellBrute + !Computes collisions in element - SUBROUTINE collision(self) + SUBROUTINE doCollisions(self) USE moduleCollisions USE moduleSpecies USE moduleList @@ -421,7 +569,9 @@ MODULE moduleMesh USE moduleRandom IMPLICIT NONE - CLASS(meshVol), INTENT(inout):: self + CLASS(meshGeneric), INTENT(inout), TARGET:: self + INTEGER:: e + CLASS(meshVol), POINTER:: vol INTEGER:: nPart !Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision INTEGER:: rnd !random index @@ -431,57 +581,57 @@ MODULE moduleMesh REAL(8):: sigmaVrelMaxNew TYPE(pointerArray), ALLOCATABLE:: partTemp(:) - nPart = self%listPart_in%amount - !Computes iterations if there is more than one particle in the cell - IF (nPart > 1) THEN - !Probability of collision - pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume + !$OMP DO SCHEDULE(DYNAMIC) + DO e=1, self%numVols + vol => self%vols(e)%obj + nPart = vol%listPart_in%amount + !Computes iterations if there is more than one particle in the cell + IF (nPart > 1) THEN + !Probability of collision + pMax = vol%totalWeight*vol%sigmaVrelMax*tauMin/vol%volume - !Increases the collisional fraction of the cell - self%collFrac = self%collFrac + REAL(nPart)*pMax*0.5D0 - - !Number of collisions in the cell - self%nColl = FLOOR(self%collFrac) + !Number of collisions in the cell + vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) - IF (self%nColl > 0) THEN - !Converts the list of particles to an array for easy access - partTemp = self%listPart_in%convert2Array() - - END IF - - DO n = 1, self%nColl - !Select random numbers - rnd = random(1, nPart) - part_i => partTemp(rnd)%part - rnd = random(1, nPart) - part_j => partTemp(rnd)%part - ij = interactionIndex(part_i%species%n, part_j%species%n) - sigmaVrelMaxNew = 0.D0 - DO k = 1, interactionMatrix(ij)%amount - CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) - - END DO - - !Update maximum cross section*v_rel per each collision - IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN - self%sigmaVrelMax = sigmaVrelMaxNew + IF (vol%nColl > 0) THEN + !Converts the list of particles to an array for easy access + partTemp = vol%listPart_in%convert2Array() END IF - !Removes one collision from the collisional fraction - self%collFrac = self%collFrac - 1.D0 - - END DO + DO n = 1, vol%nColl + !Select random numbers + rnd = random(1, nPart) + part_i => partTemp(rnd)%part + rnd = random(1, nPart) + part_j => partTemp(rnd)%part + ij = interactionIndex(part_i%species%n, part_j%species%n) + sigmaVrelMaxNew = 0.D0 + DO k = 1, interactionMatrix(ij)%amount + CALL interactionMatrix(ij)%collisions(k)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) - END IF + END DO - END SUBROUTINE collision + !Update maximum cross section*v_rel per each collision + IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN + vol%sigmaVrelMax = sigmaVrelMaxNew + + END IF + + END DO + + END IF + + END DO + !$OMP END DO + + END SUBROUTINE doCollisions !Constructs the global K matrix SUBROUTINE constructGlobalK(self) IMPLICIT NONE - CLASS(meshParticle), INTENT(inout):: self + CLASS(meshParticles), INTENT(inout):: self INTEGER:: e INTEGER, ALLOCATABLE:: n(:) REAL(8), ALLOCATABLE:: localK(:,:) diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 70bef90..73871ce 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -320,6 +320,7 @@ MODULE moduleCollisions newElectron%xi = neutral%xi newElectron%n_in = .TRUE. newElectron%vol = neutral%vol + newElectron%volColl = neutral%volColl newElectron%weight = neutral%weight newElectron%qm = electron%qm diff --git a/src/modules/moduleCompTime.f90 b/src/modules/moduleCompTime.f90 index 28273ce..d0a6755 100644 --- a/src/modules/moduleCompTime.f90 +++ b/src/modules/moduleCompTime.f90 @@ -1,11 +1,12 @@ !Information to calculate computation time MODULE moduleCompTime - REAL(8):: tStep=0.D0 - REAL(8):: tPush=0.D0 - REAL(8):: tReset=0.D0 - REAL(8):: tColl=0.D0 - REAL(8):: tWeight=0.D0 - REAL(8):: tEMField=0.D0 + REAL(8):: tStep = 0.D0 + REAL(8):: tPush = 0.D0 + REAL(8):: tReset = 0.D0 + REAL(8):: tColl = 0.D0 + REAL(8):: tCoul = 0.D0 + REAL(8):: tWeight = 0.D0 + REAL(8):: tEMField = 0.D0 END MODULE moduleCompTime diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 9660aa7..4ce168c 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -86,6 +86,7 @@ MODULE moduleInject CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) + INTEGER:: nVolColl self%id = i self%vMod = v/v_ref @@ -125,6 +126,27 @@ MODULE moduleInject IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN et = et + 1 self%edges(et) = mesh%edges(e)%obj%n + !Assign connectivity between injection edge and meshColl volume + IF (ASSOCIATED(meshForMCC, meshColl)) THEN + nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos()) + IF (nVolColl > 0) THEN + mesh%edges(e)%obj%eColl => meshColl%vols(nVolColl)%obj + + ELSE + CALL criticalError("No connection between edge and meshColl", "initInject") + + END IF + + ELSE + IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN + mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1 + + ELSE + mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2 + + END IF + + END IF END IF @@ -255,6 +277,7 @@ MODULE moduleInject CALL criticalError("No Volume associated to edge", 'addParticles') END IF + partInj(n)%volColl = randomEdge%eColl%n !Assign particle type partInj(n)%species => self%species diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 395d031..9d9de19 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -59,7 +59,7 @@ MODULE moduleInput CALL checkStatus(config, "readCase") !Read injection of particles - CALL verboseError('Reading Interactions between species...') + CALL verboseError('Reading injection of particles from boundaries...') CALL readInject(config) CALL checkStatus(config, "readInject") @@ -483,7 +483,6 @@ MODULE moduleInput END DO !Set number of particles to 0 for init state - !TODO: In a future, this should include the particles from init states nPartOld = 0 !Initialize the lock for the non-analogue (NA) list of particles @@ -497,6 +496,7 @@ MODULE moduleInput USE moduleList USE moduleCollisions USE moduleErrors + USE moduleMesh USE OMP_LIB USE json_module IMPLICIT NONE @@ -515,74 +515,98 @@ MODULE moduleInput REAL(8):: energyThreshold, energyBinding CHARACTER(:), ALLOCATABLE:: electron - CALL initInteractionMatrix(interactionMatrix) + !Firstly, checks if the object 'interactions' exists + CALL config%info('interactions', found) + IF (found) THEN + !Checks if MC collisions have been defined + CALL config%info('interactions.collisions', found) + IF (found) THEN + !Checks if a mesh for collisions has been defined + !The mesh will be initialized and reader in readGeometry + CALL config%info('interactions.meshCollisions', found) + IF (found) THEN + !Points meshForMCC to the specific mesh defined + meshForMCC => meshColl - !Path for collision cross-section data files - CALL config%get('interactions.folderCollisions', pathCollisions, found) + ELSE + !Points the meshForMCC pointer to the Particles Mesh + meshForMCC => mesh - !Inits lock for list of particles - CALL OMP_INIT_LOCK(lockCollisions) + END IF - CALL config%info('interactions.collisions', found, n_children = nInteractions) - DO i = 1, nInteractions - WRITE(iString, '(I2)') i - object = 'interactions.collisions(' // TRIM(iString) // ')' - CALL config%get(object // '.species_i', species_i, found) - pt_i = speciesName2Index(species_i) - CALL config%get(object // '.species_j', species_j, found) - pt_j = speciesName2Index(species_j) - CALL config%info(object // '.cTypes', found, n_children = nCollisions) - ij = interactionIndex(pt_i,pt_j) - !Allocates the required number of collisions per each pair of species ij - CALL interactionMatrix(ij)%init(nCollisions) + !Inits the MCC matrix + CALL initInteractionMatrix(interactionMatrix) - DO k = 1, nCollisions - WRITE (kString, '(I2)') k - object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')' - !Reads the cross section file - CALL config%get(object // '.crossSection', crossSecFile, found) - crossSecFilePath = pathCollisions // crossSecFile - IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions') - !Reads the type of collision - CALL config%get(object // '.type', cType, found) - !Initialize collision type and reads required additional data - SELECT CASE(cType) - CASE ('elastic') - !Elastic collision - CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) + !Path for collision cross-section data files + CALL config%get('interactions.folderCollisions', pathCollisions, found) - CASE ('chargeExchange') - !Resonant charge exchange - CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) + !Inits lock for list of particles + CALL OMP_INIT_LOCK(lockCollisions) - CASE ('ionization') - !Electorn impact ionization - CALL config%get(object // '.energyThreshold', energyThreshold, found) - IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions') - CALL config%get(object // '.electron', electron, found) - IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') - CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron) + CALL config%info('interactions.collisions', found, n_children = nInteractions) + DO i = 1, nInteractions + WRITE(iString, '(I2)') i + object = 'interactions.collisions(' // TRIM(iString) // ')' + CALL config%get(object // '.species_i', species_i, found) + pt_i = speciesName2Index(species_i) + CALL config%get(object // '.species_j', species_j, found) + pt_j = speciesName2Index(species_j) + CALL config%info(object // '.cTypes', found, n_children = nCollisions) + ij = interactionIndex(pt_i,pt_j) + !Allocates the required number of collisions per each pair of species ij + CALL interactionMatrix(ij)%init(nCollisions) - CASE ('recombination') - !Electorn impact ionization - CALL config%get(object // '.energyBinding', energyBinding, found) - IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions') - CALL config%get(object // '.electron', electron, found) - IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') - CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, & - crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron) + DO k = 1, nCollisions + WRITE (kString, '(I2)') k + object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')' + !Reads the cross section file + CALL config%get(object // '.crossSection', crossSecFile, found) + crossSecFilePath = pathCollisions // crossSecFile + IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions') + !Reads the type of collision + CALL config%get(object // '.type', cType, found) + !Initialize collision type and reads required additional data + SELECT CASE(cType) + CASE ('elastic') + !Elastic collision + CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, & + crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) - CASE DEFAULT - CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions') + CASE ('chargeExchange') + !Resonant charge exchange + CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, & + crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) - END SELECT + CASE ('ionization') + !Electorn impact ionization + CALL config%get(object // '.energyThreshold', energyThreshold, found) + IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions') + CALL config%get(object // '.electron', electron, found) + IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') + CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, & + crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron) - END DO + CASE ('recombination') + !Electorn impact ionization + CALL config%get(object // '.energyBinding', energyBinding, found) + IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions') + CALL config%get(object // '.electron', electron, found) + IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions') + CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, & + crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron) - END DO + CASE DEFAULT + CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions') + + END SELECT + + END DO + + END DO + + END IF + + END IF END SUBROUTINE readInteractions @@ -702,17 +726,23 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config LOGICAL:: found + LOGICAL:: doubleMesh CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile CHARACTER(:), ALLOCATABLE:: fullPath + !Firstly, indicates if a specific mesh for MC collisions is being use + doubleMesh = ASSOCIATED(meshForMCC, meshColl) + !Selects the type of geometry. CALL config%get('geometry.type', mesh%geometry, found) + IF (doubleMesh) meshColl%geometry = mesh%geometry !Gets the type of mesh CALL config%get('geometry.meshType', meshFormat, found) SELECT CASE(meshFormat) CASE ("gmsh2") CALL initGmsh2(mesh) + IF (doubleMesh) CALL initGmsh2(meshColl) CASE DEFAULT CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry") @@ -723,6 +753,14 @@ MODULE moduleInput CALL config%get('geometry.meshFile', meshFile, found) fullpath = path // meshFile CALL mesh%readMesh(fullPath) + DEALLOCATE(fullPath, meshFile) + IF (doubleMesh) THEN + !Reads the mesh file for collisions + CALL config%get('interactions.meshCollisions', meshFile, found) + fullpath = path // meshFile + CALL meshColl%readMesh(fullPath) + + END IF !Creates the connectivity between elements SELECT CASE(mesh%geometry) @@ -744,9 +782,24 @@ MODULE moduleInput END SELECT CALL mesh%connectMesh - !Builds the K matrix + IF (doubleMesh) THEN + meshColl%connectMesh => mesh%connectMesh + CALL meshColl%connectMesh + + END IF + + !Builds the K matrix for the Particles mesh CALL mesh%constructGlobalK() + !Assign the procedure to find a volume for meshColl + IF (doubleMesh) THEN + findCellColl => findCellCollMesh + + ELSE + findCellColl => findCellSameMesh + + END IF + END SUBROUTINE readGeometry SUBROUTINE readEMBoundary(config) diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index c08061a..acdd1de 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -10,7 +10,7 @@ MODULE moduleList END TYPE lNode TYPE listNode - INTEGER:: amount = 0!TODO: Make private + INTEGER:: amount = 0 TYPE(lNode),POINTER:: head => NULL() TYPE(lNode),POINTER:: tail => NULL() CONTAINS diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90 index 963302d..a9efdc4 100644 --- a/src/modules/moduleOutput.f90 +++ b/src/modules/moduleOutput.f90 @@ -109,7 +109,7 @@ MODULE moduleOutput END IF - WRITE (20, "(I10, I10, 6(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight, tEMField + WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField CLOSE(20) diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index b549104..089ae71 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -483,21 +483,6 @@ MODULE moduleSolver END SUBROUTINE push1DRadCharged - !Do the collisions in all the cells - SUBROUTINE doCollisions() - USE moduleMesh - IMPLICIT NONE - - INTEGER:: e - - !$OMP DO SCHEDULE(DYNAMIC) - DO e=1, mesh%numVols - CALL mesh%vols(e)%obj%collision() - END DO - !$OMP END DO - - END SUBROUTINE doCollisions - SUBROUTINE doReset() USE moduleSpecies USE moduleMesh @@ -623,6 +608,14 @@ MODULE moduleSolver END DO + !$OMP SECTION + !Erase the list of particles inside the cell in coll mesh + DO e = 1, meshColl%numVols + meshColl%vols(e)%obj%totalWeight = 0.D0 + CALL meshColl%vols(e)%obj%listPart_in%erase() + + END DO + !$OMP END SECTIONS !$OMP SINGLE @@ -788,6 +781,7 @@ MODULE moduleSolver volOld => mesh%vols(part%vol)%obj CALL volOld%findCell(part) + CALL findCellColl(part) volNew => mesh%vols(part%vol)%obj !Call the NA shcme IF (ASSOCIATED(self%weightingScheme)) THEN @@ -831,7 +825,7 @@ MODULE moduleSolver counterOutput=0 CALL mesh%printOutput(t) - CALL mesh%printColl(t) + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) CALL mesh%printEM(t) WRITE(*, "(5X,A21,I10,A1,I10)") "t/tmax: ", t, "/", tmax WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index 7d89e3b..39c1b20 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -39,6 +39,7 @@ MODULE moduleSpecies REAL(8):: v(1:3) !Velocity CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle INTEGER:: vol !Index of element in which the particle is located + INTEGER:: volColl !Index of element in which the particle is located in the Collision Mesh 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