From 73fc9f69c1d22022e99dceebb63dd6a63ad38c82 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 25 Oct 2020 08:08:18 +0100 Subject: [PATCH] Main changes: - Injection is now performed in parallalel (an IF statement could be required to avoid overhead when number of injected particles is below a margin). - Added the possibility for multiple injections. --- src/DSMC_Neutrals.f95 | 3 +- src/modules/moduleInject.f95 | 53 ++++++++--------------------------- src/modules/moduleList.f95 | 6 ++++ src/modules/moduleMesh.f95 | 19 +++++++++++++ src/modules/moduleMeshCyl.f95 | 30 +++++++++++--------- src/modules/moduleSolver.f95 | 45 ++++++++++++++++------------- 6 files changed, 80 insertions(+), 76 deletions(-) diff --git a/src/DSMC_Neutrals.f95 b/src/DSMC_Neutrals.f95 index 50f4f29..f71bdbb 100644 --- a/src/DSMC_Neutrals.f95 +++ b/src/DSMC_Neutrals.f95 @@ -37,11 +37,10 @@ PROGRAM DSMC_Neutrals !$OMP SINGLE tStep = omp_get_wtime() tPush = omp_get_wtime() + !$OMP END SINGLE CALL doInjects() - !$OMP END SINGLE NOWAIT - !Push old particles CALL doPushes() diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95 index 63f53fb..9b87133 100644 --- a/src/modules/moduleInject.f95 +++ b/src/modules/moduleInject.f95 @@ -119,23 +119,19 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(in):: self REAL(8):: randomX INTEGER:: j - INTEGER:: randomEdge REAL(8):: randomPos - !Edge nodes - INTEGER:: n1 = 0, n2 = 0 - !Edge nodes coordinates - REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0 INTEGER:: nMin, nMax !Min and Max index in partInj array INTEGER:: n + CLASS(meshEdge), POINTER:: randomEdge !Insert particles - !TODO: Adjust for multiple injections - nMin = 1 - nMax = self%nParticles + nMin = SUM(inject(1:(self%id-1))%nParticles) + 1 + nMax = nMin + self%nParticles !Assign particle type partInj(nMin:nMax)%pt = self%pt !Assign weight to particle. partInj(nMin:nMax)%weight = species(self%pt)%obj%weight + !$OMP DO DO n = nMin, nMax !Select edge randomly from which inject particle randomX = RAND()*self%sumWeight @@ -145,34 +141,16 @@ MODULE moduleInject END DO - randomEdge = self%edges(j) - - !Get coordinates of edge nodes - SELECT TYPE(edge => mesh%edges(randomEdge)%obj) - CLASS IS (meshEdgeCyl) - n1 = edge%n1%n - n2 = edge%n2%n - - END SELECT - SELECT TYPE(node => mesh%nodes(n1)%obj) - TYPE IS (meshNodeCyl) - p1(1) = node%z - p1(2) = node%r - - END SELECT - SELECT TYPE(node => mesh%nodes(n2)%obj) - TYPE IS (meshNodeCyl) - p2(1) = node%z - p2(2) = node%r - - END SELECT - + randomEdge => mesh%edges(self%edges(j))%obj + !Random position in edge + randomPos = RAND() + partInj(n)%r = randomEdge%randPos(randomPos) !Volume associated to the edge: - IF (DOT_PRODUCT(self%n, mesh%edges(randomEdge)%obj%normal) >= 0.D0) THEN - partInj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n + IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN + partInj(n)%e_p = randomEdge%e1%n ELSE - partInj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n + partInj(n)%e_p = randomEdge%e2%n END IF @@ -180,20 +158,13 @@ MODULE moduleInject vBC(self%v(2), self%vTh(2)), & vBC(self%v(3), self%vTh(3)) /) - !Random position in edge - !TODO: Use edge coordinates and transformations for this process - randomPos = RAND() - partInj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1)) - partInj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2)) - partInj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3)) - - !Push new particle CALL push(partInj(n)) !Assign cell to new particle CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n)) END DO + !$OMP END DO END SUBROUTINE addParticlesMaxwellian diff --git a/src/modules/moduleList.f95 b/src/modules/moduleList.f95 index 0abdd5a..5dabe30 100644 --- a/src/modules/moduleList.f95 +++ b/src/modules/moduleList.f95 @@ -1,6 +1,7 @@ !Linked list of particles MODULE moduleList USE moduleSpecies + USE OMP_LIB IMPLICIT NONE TYPE lNode @@ -20,6 +21,10 @@ MODULE moduleList END TYPE listNode + INTEGER(KIND=OMP_LOCK_KIND):: resetLock + + TYPE(listNode):: partNew + TYPE pointerArray TYPE(particle), POINTER:: part @@ -49,6 +54,7 @@ MODULE moduleList END SUBROUTINE addToList + !converts list to array FUNCTION convert2Array(self) RESULT(partArray) IMPLICIT NONE diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95 index df637f2..85048ba 100644 --- a/src/modules/moduleMesh.f95 +++ b/src/modules/moduleMesh.f95 @@ -59,6 +59,8 @@ MODULE moduleMesh INTEGER:: bt = 0 CONTAINS PROCEDURE(initEdge_interface), DEFERRED, PASS:: init + PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary + PROCEDURE(randPos_interface), DEFERRED, PASS:: randPos END TYPE meshEdge @@ -73,6 +75,23 @@ MODULE moduleMesh END SUBROUTINE initEdge_interface + SUBROUTINE boundary_interface(self, part) + USE moduleSpecies + + IMPORT:: meshEdge + CLASS (meshEdge), INTENT(inout):: self + CLASS (particle), INTENT(inout):: part + + END SUBROUTINE + + PURE FUNCTION randPos_interface(self, rnd) RESULT(r) + IMPORT:: meshEdge + CLASS(meshEdge), INTENT(in):: self + REAL(8), INTENT(in):: rnd + REAL(8):: r(1:3) + + END FUNCTION randPos_interface + END INTERFACE !Containers for edges in the mesh diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95 index 2175119..69c9bc6 100644 --- a/src/modules/moduleMeshCyl.f95 +++ b/src/modules/moduleMeshCyl.f95 @@ -25,23 +25,11 @@ MODULE moduleMeshCyl !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS - PROCEDURE, PASS:: init => initEdgeCyl - PROCEDURE(boundary_interface), PASS, DEFERRED:: fBoundary + PROCEDURE, PASS:: init => initEdgeCyl + PROCEDURE, PASS:: randPos => randPosCyl END TYPE meshEdgeCyl - ABSTRACT INTERFACE - SUBROUTINE boundary_interface(self, part) - USE moduleSpecies - - IMPORT:: meshEdgeCyl - CLASS (meshEdgeCyl), INTENT(inout):: self - CLASS (particle), INTENT(inout):: part - - END SUBROUTINE - - END INTERFACE - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl CONTAINS PROCEDURE, PASS:: collision => collision2DCyl @@ -140,6 +128,20 @@ MODULE moduleMeshCyl END SUBROUTINE initEdgeCyl + + PURE FUNCTION randPosCyl(self, rnd) RESULT(r) + CLASS(meshEdgeCyl), INTENT(in):: self + REAL(8), INTENT(in):: rnd + REAL(8):: r(1:3) + REAL(8):: p1(1:2), p2(1:2) + + p1 = (/self%z(1), self%r(1) /) + p2 = (/self%z(2), self%r(2) /) + r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 + r(3) = 0.D0 + + END FUNCTION randPosCyl + !Vol functions !Quad Volume SUBROUTINE initVolQuadCyl(self, n, p) diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95 index 1676425..ec05cd7 100644 --- a/src/modules/moduleSolver.f95 +++ b/src/modules/moduleSolver.f95 @@ -74,31 +74,14 @@ MODULE moduleSolver END SUBROUTINE doCollisions - !Scatter particles in the grid - SUBROUTINE doScatter - USE moduleSpecies - USE moduleMesh - - INTEGER:: n - - !Loops over the particles to scatter them - !$OMP DO - DO n=1, nPartOld - CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n)) - - END DO - !$OMP END DO - - END SUBROUTINE doScatter - SUBROUTINE doReset() USE moduleSpecies IMPLICIT NONE INTEGER:: nn, n - INTEGER:: nPartNew + INTEGER, SAVE:: nPartNew INTEGER, SAVE:: nInjIn, nOldIn - TYPE(particle), ALLOCATABLE:: partTemp(:) + TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) !$OMP SECTIONS !$OMP SECTION @@ -121,7 +104,10 @@ MODULE moduleSolver CALL MOVE_ALLOC(partOld, partTemp) nPartNew = nInjIn + nOldIn ALLOCATE(partOld(1:nPartNew)) + !$OMP END SINGLE + !$OMP SECTIONS + !$OMP SECTION nn = 0 DO n = 1, nPartInj IF (partInj(n)%n_in) THEN @@ -132,6 +118,8 @@ MODULE moduleSolver END DO + !$OMP SECTION + nn = nInjIn DO n = 1, nPartOld IF (partTemp(n)%n_in) THEN nn = nn + 1 @@ -140,12 +128,31 @@ MODULE moduleSolver END IF END DO + !$OMP END SECTIONS + !$OMP SINGLE nPartOld = nPartNew !$OMP END SINGLE END SUBROUTINE doReset + !Scatter particles in the grid + SUBROUTINE doScatter + USE moduleSpecies + USE moduleMesh + + INTEGER:: n + + !Loops over the particles to scatter them + !$OMP DO + DO n=1, nPartOld + CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n)) + + END DO + !$OMP END DO + + END SUBROUTINE doScatter + SUBROUTINE doOutput(t) USE moduleMesh USE moduleOutput