diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 0237de1..dacbf9c 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -40,10 +40,10 @@ MODULE moduleRandom INTEGER:: rnd REAL(8):: rnd01 - rnd = 0.D0 + rnd = 0 CALL RANDOM_NUMBER(rnd01) - rnd = INT(REAL(b - a) * rnd01) + 1 + rnd = NINT(REAL(b - a) * rnd01) + a END FUNCTION randomIntAB diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 90773e5..2d2d6c6 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -235,6 +235,7 @@ MODULE moduleMesh2DCyl p1 = (/self%z(1), self%r(1) /) p2 = (/self%z(2), self%r(2) /) r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 + r(2) = (self%r(2) + self%r(1)) * 0.5D0 r(3) = 0.D0 END FUNCTION randPosEdge @@ -578,9 +579,6 @@ MODULE moduleMesh2DCyl REAL(8):: detJ REAL(8):: fPsi(1:4) REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) - REAL(8):: nodeR(1:3) - REAL(8), PARAMETER:: correction = 0.81D0 !TODO: No idea why this is needed, this needs more attention. Probably a better - ! calculation of the volume is needed. self%volume = 0.D0 @@ -594,26 +592,18 @@ MODULE moduleMesh2DCyl !Computes total volume of the cell self%volume = r*detJ*PI8 !4*2*pi !Computes volume per node - self%n1%v = self%n1%v + fPsi(1)*self%volume - nodeR = self%n1%getCoordinates() - if (nodeR(2) == 0.D0) then - self%n1%v = self%n1%v * correction - end if - self%n2%v = self%n2%v + fPsi(2)*self%volume - nodeR = self%n2%getCoordinates() - if (nodeR(2) == 0.D0) then - self%n2%v = self%n2%v * correction - end if - self%n3%v = self%n3%v + fPsi(3)*self%volume - nodeR = self%n3%getCoordinates() - if (nodeR(2) == 0.D0) then - self%n3%v = self%n3%v * correction - end if - self%n4%v = self%n4%v + fPsi(4)*self%volume - nodeR = self%n4%getCoordinates() - if (nodeR(2) == 0.D0) then - self%n4%v = self%n4%v * correction - end if + 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 END SUBROUTINE volumeQuad diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index fb01f75..57e3131 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -283,6 +283,7 @@ MODULE moduleInject USE moduleMesh USE moduleRandom USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self @@ -312,12 +313,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 * sqrt(partInj(n)%r(2) / (0.1D0*L_ref)) !Volume associated to the edge: IF (ASSOCIATED(randomEdge%e1)) THEN partInj(n)%cell = randomEdge%e1%n