Merge branch 'issue/injection' into 'development'

Change in injection to achieve uniform density of particles.

See merge request JorgeGonz/fpakc!50

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
 

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
 

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