From 6389c8ba2d090db7cc4b29f4c775f4c18f8395d9 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 27 Jun 2024 12:08:08 +0200 Subject: [PATCH 1/6] Quick because food Cartesian fixed now --- src/modules/common/moduleRandom.f90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index cd553a8..3fbf5c8 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -73,10 +73,21 @@ MODULE moduleRandom REAL(8), INTENT(in):: cumWeight(1:) REAL(8), INTENT(in):: sumWeight REAL(8):: rnd0b - INTEGER:: rnd + INTEGER:: rnd, i rnd0b = random(0.D0, sumWeight) - rnd = MINLOC(DABS(rnd0b - cumWeight), 1) + i = 1 + DO + IF (rnd0b <= cumWeight(i)) THEN + rnd = i + EXIT + + ELSE + i = i +1 + + END IF + END DO + ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1) END FUNCTION randomWeighted -- 2.49.1 From 5386114d15114574807f88cac4808924537e3c10 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 29 Jun 2024 14:58:48 +0200 Subject: [PATCH 2/6] Cylindrical injection working better Seems things are a bit better. Still, more cases are needed and still not perfectly uniform... --- src/modules/common/moduleRandom.f90 | 4 +-- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 37 +++++++++++++--------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 3fbf5c8..0237de1 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -75,10 +75,10 @@ MODULE moduleRandom REAL(8):: rnd0b INTEGER:: rnd, i - rnd0b = random(0.D0, sumWeight) + rnd0b = random() i = 1 DO - IF (rnd0b <= cumWeight(i)) THEN + IF (rnd0b <= cumWeight(i)/sumWeight) THEN rnd = i EXIT diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index dfa3035..dbe182f 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = r2(2)**2 - r1(2)**2 + self%weight = DABS(self%r(2)**2 - self%r(1)**2) !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & @@ -223,21 +223,12 @@ MODULE moduleMesh2DCyl CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd REAL(8):: r(1:3) - REAL(8):: dr, dz + REAL(8):: p1(1:2), p2(1:2) rnd = random() - dr = self%r(2) - self%r(1) - dz = self%z(2) - self%z(1) - IF (dr /= 0.D0) THEN - r(2) = dr * DSQRT(rnd) + self%r(1) - r(1) = dz * (r(2) - self%r(1))/dr + self%z(1) - - ELSE - r(2) = self%r(1) - r(1) = dz * rnd + self%z(1) - - END IF - + 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 randPosEdge @@ -572,7 +563,7 @@ MODULE moduleMesh2DCyl !Compute element volume PURE SUBROUTINE volumeQuad(self) - USE moduleConstParam, ONLY: PI8 + USE moduleConstParam, ONLY: PI, PI8 IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(inout):: self @@ -596,15 +587,31 @@ MODULE moduleMesh2DCyl 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 + IF (self%r(1) == 0.D0) THEN + self%n1%v = self%n1%v * PI / 2.D0 + + END IF 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 + IF (self%r(2) == 0.D0) THEN + self%n2%v = self%n2%v * PI / 2.D0 + + END IF 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 + IF (self%r(3) == 0.D0) THEN + self%n3%v = self%n3%v * PI / 2.D0 + + END IF 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 + IF (self%r(4) == 0.D0) THEN + self%n4%v = self%n4%v * PI / 2.D0 + + END IF END SUBROUTINE volumeQuad -- 2.49.1 From 4cadfe536781bd7b901e7adaf8a361be7ebc9c13 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 29 Jun 2024 22:22:10 +0200 Subject: [PATCH 3/6] Seems are a bit better There is still less density in the axis. I don't find a reason why. There must be a modification to the weight... --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 24 ++++------------------ src/modules/moduleInject.f90 | 5 ++--- 2 files changed, 6 insertions(+), 23 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index dbe182f..4991b7e 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -586,32 +586,16 @@ MODULE moduleMesh2DCyl !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 - IF (self%r(1) == 0.D0) THEN - self%n1%v = self%n1%v * PI / 2.D0 - - END IF + self%n1%v = self%n1%v + fPsi(1)*self%volume 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 - IF (self%r(2) == 0.D0) THEN - self%n2%v = self%n2%v * PI / 2.D0 - - END IF + self%n2%v = self%n2%v + fPsi(2)*self%volume 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 - IF (self%r(3) == 0.D0) THEN - self%n3%v = self%n3%v * PI / 2.D0 - - END IF + self%n3%v = self%n3%v + fPsi(3)*self%volume 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 - IF (self%r(4) == 0.D0) THEN - self%n4%v = self%n4%v * PI / 2.D0 - - END IF + self%n4%v = self%n4%v + fPsi(4)*self%volume END SUBROUTINE volumeQuad diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index b37e021..fb01f75 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -306,8 +306,6 @@ MODULE moduleInject nMin = nMin + 1 nMax = nMin + self%nParticles - 1 - !Assign weight to particle. - partInj(nMin:nMax)%weight = self%species%weight !Particle is considered to be outside the domain partInj(nMin:nMax)%n_in = .FALSE. !$OMP END SINGLE @@ -315,10 +313,11 @@ MODULE moduleInject !$OMP DO DO n = nMin, nMax randomX = randomWeighted(self%cumWeight, self%sumWeight) - 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 !Volume associated to the edge: IF (ASSOCIATED(randomEdge%e1)) THEN partInj(n)%cell = randomEdge%e1%n -- 2.49.1 From cd7bf66bd84213b6cd385dbe84fb7edd6cd0a5eb Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 30 Jun 2024 10:36:36 +0200 Subject: [PATCH 4/6] A workaround The random position for edges in the axis is corrected so that there is a more uniform charge density in the axis. Still, things are not perfect and this is something to really look into in the future. --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 4991b7e..9278a27 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -226,6 +226,12 @@ MODULE moduleMesh2DCyl REAL(8):: p1(1:2), p2(1:2) rnd = random() + IF (self%r(1) == 0.D0 .OR. & + self%r(2) == 0.D0) THEN + rnd = rnd**(1.D0/3.D0) + + END IF + p1 = (/self%z(1), self%r(1) /) p2 = (/self%z(2), self%r(2) /) r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 -- 2.49.1 From 59a322a4c7df17c6de8d507aef35588eae28935d Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 30 Jun 2024 10:46:05 +0200 Subject: [PATCH 5/6] Clean up Fixing calculation of node volumes. --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 9278a27..df01782 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -580,6 +580,7 @@ MODULE moduleMesh2DCyl REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) self%volume = 0.D0 + !2D 1 point Gauss Quad Integral Xi = 0.D0 dPsi = self%dPsi(Xi, 4) @@ -590,17 +591,9 @@ 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)*self%volume - Xi = (/ 5.D-1, -5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) self%n2%v = self%n2%v + fPsi(2)*self%volume - Xi = (/ 5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) self%n3%v = self%n3%v + fPsi(3)*self%volume - Xi = (/-5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) self%n4%v = self%n4%v + fPsi(4)*self%volume END SUBROUTINE volumeQuad -- 2.49.1 From e277fe6ddb8526ec2c1ae8fedccd4a777c178c0d Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 30 Jun 2024 08:48:51 +0000 Subject: [PATCH 6/6] PI is no needed here --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index df01782..7bae6b3 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -569,7 +569,7 @@ MODULE moduleMesh2DCyl !Compute element volume PURE SUBROUTINE volumeQuad(self) - USE moduleConstParam, ONLY: PI, PI8 + USE moduleConstParam, ONLY: PI8 IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(inout):: self -- 2.49.1