From b36f9c26152339aef154402ebabce74f4075d866 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 9 Jul 2024 17:49:42 +0200 Subject: [PATCH] Shifting towards constant number of particles per edge So now each edge has the same number of particles and the weight of each particle is calculated based on the surface of each edge compared to the total one. Only in 2DCyl, still to extend to other geometries. Not perfect constant density, but the issue might be the node volume. --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 27 ++++++++++--------- .../mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 | 14 ---------- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 26 +++++++++--------- src/modules/mesh/moduleMesh.f90 | 2 ++ src/modules/moduleInject.f90 | 16 +++++------ 5 files changed, 36 insertions(+), 49 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 426846d..3f88902 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -144,6 +144,7 @@ MODULE moduleMesh2DCyl USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleConstParam, ONLY: PI IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(out):: self @@ -163,7 +164,15 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = DABS(self%r(2)**2 - self%r(1)**2) + !Edge surface + IF (self%z(2) /= self%z(1)) THEN + self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1)) + + ELSE + self%surface = ABS(self%r(2)**2 - self%r(1)**2) + + END IF + self%surface = self%surface * PI !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & @@ -586,18 +595,10 @@ MODULE moduleMesh2DCyl !Computes total volume of the cell self%volume = r*detJ*PI8 !4*2*pi !Computes volume per node - Xi = (/-5.D-1, -5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) - self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 - Xi = (/ 5.D-1, -5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) - self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 - Xi = (/ 5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) - self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 - Xi = (/-5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) - self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 + self%n1%v = self%n1%v + fPsi(1)*self%volume + self%n2%v = self%n2%v + fPsi(2)*self%volume + self%n3%v = self%n3%v + fPsi(3)*self%volume + self%n4%v = self%n4%v + fPsi(4)*self%volume END SUBROUTINE volumeQuad diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 2c02358..0df1289 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -296,20 +296,6 @@ MODULE moduleMeshInputGmsh2 CLOSE(10) - ! Adjust node volume at axis - SELECT CASE(self%geometry) - CASE("Cyl") - DO n = 1, self%numNodes - r = self%nodes(n)%obj%getCoordinates() - IF (r(2) == 0.D0) THEN - self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0 - - END IF - - END DO - - END SELECT - !Call mesh connectivity CALL self%connectMesh diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index eb6af79..309ca5c 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -496,19 +496,19 @@ MODULE moduleMeshInputVTU END DO - ! Adjust node volume at axis - SELECT CASE(self%geometry) - CASE("Cyl") - DO n = 1, numNodes - r = self%nodes(n)%obj%getCoordinates() - IF (r(2) == 0.D0) THEN - self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0 - - END IF - - END DO - - END SELECT + ! ! Adjust node volume at axis + ! SELECT CASE(self%geometry) + ! CASE("Cyl") + ! DO n = 1, numNodes + ! r = self%nodes(n)%obj%getCoordinates() + ! IF (r(2) == 0.D0) THEN + ! self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0 + ! + ! END IF + ! + ! END DO + ! + ! END SELECT !Call mesh connectivity CALL self%connectMesh diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index e96ff2a..af1cd36 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -78,6 +78,8 @@ MODULE moduleMesh REAL(8):: normal(1:3) !Weight for random injection of particles REAL(8):: weight = 1.D0 + ! Surface of edge + REAL(8):: surface = 0.D0 !Pointer to boundary type TYPE(boundaryCont), POINTER:: boundary !Array of functions for boundary conditions diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 69a3dd6..1d4e2e7 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -63,6 +63,7 @@ MODULE moduleInject INTEGER, ALLOCATABLE:: edges(:) !Array with edges REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability REAL(8):: sumWeight + REAL(8):: surface ! Total surface of injection TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction CONTAINS PROCEDURE, PASS:: init => initInject @@ -164,15 +165,12 @@ MODULE moduleInject END DO - !Calculates cumulative probability - ALLOCATE(self%cumWeight(1:self%nEdges)) - et = 1 - self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight - DO et = 2, self%nEdges - self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1) + !Calculates total area + self%surface = 0.D0 + DO et = 1, self%nEdges + self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface END DO - self%sumWeight = self%cumWeight(self%nEdges) END SUBROUTINE initInject @@ -313,12 +311,12 @@ MODULE moduleInject !$OMP DO DO n = nMin, nMax - randomX = randomWeighted(self%cumWeight, self%sumWeight) + randomX = random(1, self%nEdges) randomEdge => mesh%edges(self%edges(randomX))%obj !Random position in edge partInj(n)%r = randomEdge%randPos() !Assign weight to particle. - partInj(n)%weight = self%species%weight + partInj(n)%weight = self%species%weight * randomEdge%surface / self%surface !Volume associated to the edge: IF (ASSOCIATED(randomEdge%e1)) THEN partInj(n)%cell = randomEdge%e1%n