diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index cd553a8..0237de1 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) + rnd0b = random() + i = 1 + DO + IF (rnd0b <= cumWeight(i)/sumWeight) THEN + rnd = i + EXIT + + ELSE + i = i +1 + + END IF + END DO + ! rnd = MINLOC(DABS(rnd0b - cumWeight), 1) END FUNCTION randomWeighted diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index dfa3035..7bae6b3 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,18 @@ 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) + 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 r(3) = 0.D0 END FUNCTION randPosEdge @@ -583,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) @@ -593,18 +591,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/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