First commit of branch performance:
Bugs fixed: - Solved an issue with particles being injected with infinite velocity resulting in Inf velocity in some cells of the output files. - Particles are now equally distributed in cylindrical geometry along the radial direction. New features: - Particles now have their own weight that is recalculated when the particle moves to a new cell. This avoid the reduction of density at r = 0. Cases: - Added a case of Argon flow around a cylinder to measure performance and future improvements.
This commit is contained in:
parent
bd7e8b040b
commit
05f5adcfe1
10 changed files with 5003 additions and 28 deletions
|
|
@ -9,7 +9,10 @@ MODULE moduleInject
|
|||
REAL(8):: n(1:3) !Direction of injection
|
||||
INTEGER:: nParticles !Number of particles to introduce each time step
|
||||
INTEGER:: pt !Species of injection
|
||||
INTEGER:: nEdges
|
||||
INTEGER, ALLOCATABLE:: edges(:) !Array with edges
|
||||
REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection
|
||||
REAL(8):: sumWeight
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: init => initInject
|
||||
PROCEDURE, PASS:: addParticles => addParticlesMaxwellian
|
||||
|
|
@ -22,6 +25,7 @@ MODULE moduleInject
|
|||
CONTAINS
|
||||
SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface)
|
||||
USE moduleMesh
|
||||
USE moduleMeshCyl
|
||||
USE moduleRefParam
|
||||
USE moduleConstParam
|
||||
USE moduleSpecies
|
||||
|
|
@ -32,7 +36,7 @@ MODULE moduleInject
|
|||
REAL(8), INTENT(in):: v, n(1:3), T(1:3)
|
||||
INTEGER, INTENT(in):: pt, physicalSurface
|
||||
REAL(8), INTENT(in):: flow
|
||||
INTEGER:: nEdges, e, et
|
||||
INTEGER:: e, et
|
||||
INTEGER:: phSurface(1:mesh%numEdges)
|
||||
|
||||
self%id = i
|
||||
|
|
@ -49,17 +53,24 @@ MODULE moduleInject
|
|||
|
||||
END DO
|
||||
|
||||
nEdges = COUNT(phSurface == physicalSurface)
|
||||
ALLOCATE(inject(i)%edges(1:nEdges))
|
||||
self%nEdges = COUNT(phSurface == physicalSurface)
|
||||
ALLOCATE(inject(i)%edges(1:self%nEdges))
|
||||
ALLOCATE(inject(i)%weight(1:self%nEdges))
|
||||
et = 0
|
||||
DO e=1, mesh%numEdges
|
||||
IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
|
||||
et = et + 1
|
||||
self%edges(et) = mesh%edges(e)%obj%n
|
||||
self%edges(et) = mesh%edges(e)%obj%n
|
||||
SELECT TYPE(edge => mesh%edges(e)%obj)
|
||||
CLASS IS (meshEdgeCyl)
|
||||
self%weight(et) = (edge%r(1)+edge%r(2))/2.D0
|
||||
|
||||
END SELECT
|
||||
|
||||
END IF
|
||||
|
||||
END DO
|
||||
self%sumWeight = SUM(self%weight)
|
||||
nPartInj = nPartInj + self%nParticles
|
||||
|
||||
END SUBROUTINE
|
||||
|
|
@ -70,8 +81,11 @@ MODULE moduleInject
|
|||
REAL(8), INTENT (in):: u, vth
|
||||
REAL(8):: x, y
|
||||
vBC=0.D0
|
||||
x = 0.D0
|
||||
|
||||
x=RAND()
|
||||
DO WHILE (x == 0.D0)
|
||||
x=RAND()
|
||||
END DO
|
||||
y=RAND()
|
||||
|
||||
vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
|
||||
|
|
@ -86,6 +100,8 @@ MODULE moduleInject
|
|||
IMPLICIT NONE
|
||||
|
||||
CLASS(injectGeneric), INTENT(in):: self
|
||||
REAL(8):: randomX
|
||||
INTEGER:: j
|
||||
INTEGER:: randomEdge
|
||||
REAL(8):: randomPos
|
||||
REAL(8):: vVec(1:3), vTh(1:3)
|
||||
|
|
@ -93,15 +109,31 @@ MODULE moduleInject
|
|||
INTEGER:: n1 = 0, n2 = 0
|
||||
!Edge nodes coordinates
|
||||
REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0
|
||||
INTEGER:: nMin, nMax !Min and Max index in part_inj array
|
||||
INTEGER:: n
|
||||
|
||||
vVec = self%v * self%n
|
||||
vTh = DSQRT(self%T/species(self%pt)%obj%m)
|
||||
|
||||
!Insert particles
|
||||
DO n = 1, self%nParticles
|
||||
!TODO: Adjust for multiple injections
|
||||
nMin = 1
|
||||
nMax = self%nParticles
|
||||
!Assign particle type
|
||||
part_inj(nMin:nMax)%pt = self%pt
|
||||
!Assign weight to particle.
|
||||
part_inj(nMin:nMax)%weight = species(self%pt)%obj%weight
|
||||
part_inj(nMin:nMax)%n_in = .TRUE.
|
||||
DO n = nMin, nMax
|
||||
!Select edge randomly from which inject particle
|
||||
randomEdge = self%edges(INT(RAND()*(SIZE(self%edges)-1)+1))
|
||||
randomX = RAND()*self%sumWeight
|
||||
DO j = 1, self%nEdges
|
||||
IF (randomX < self%weight(j)) EXIT
|
||||
randomX = randomX - self%weight(j)
|
||||
|
||||
END DO
|
||||
|
||||
randomEdge = self%edges(j)
|
||||
|
||||
!Get coordinates of edge nodes
|
||||
SELECT TYPE(edge => mesh%edges(randomEdge)%obj)
|
||||
|
|
@ -132,8 +164,6 @@ MODULE moduleInject
|
|||
|
||||
END IF
|
||||
|
||||
part_inj(n)%pt = self%pt
|
||||
part_inj(n)%n_in = .TRUE.
|
||||
part_inj(n)%v = (/ vBC(vVec(1), vTh(1)), &
|
||||
vBC(vVec(2), vTh(2)), &
|
||||
vBC(vVec(3), vTh(3)) /)
|
||||
|
|
@ -145,6 +175,7 @@ MODULE moduleInject
|
|||
part_inj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2))
|
||||
part_inj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3))
|
||||
|
||||
|
||||
!Push new particle
|
||||
CALL push(part_inj(n))
|
||||
|
||||
|
|
|
|||
|
|
@ -74,11 +74,11 @@ MODULE moduleInput
|
|||
IF (.NOT. found) CALL criticalError('Reference radius not found','readReference')
|
||||
|
||||
!Derived parameters
|
||||
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
|
||||
sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section
|
||||
L_ref = 1.D0/(sigma_ref*n_ref) !mean free path
|
||||
ti_ref = L_ref/v_ref !reference time
|
||||
Vol_ref = L_ref**3 !reference volume
|
||||
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
|
||||
|
||||
END SUBROUTINE readReference
|
||||
|
||||
|
|
@ -104,7 +104,10 @@ MODULE moduleInput
|
|||
IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
|
||||
|
||||
!Convert simulation time to number of iterations
|
||||
tmax = INT(time/(ti_ref*tau))
|
||||
tmax = INT(time/tau)
|
||||
|
||||
!Makes tau non-dimensional
|
||||
tau = tau / ti_ref
|
||||
|
||||
END SUBROUTINE readCase
|
||||
|
||||
|
|
|
|||
|
|
@ -124,11 +124,12 @@ MODULE moduleMesh
|
|||
|
||||
END SUBROUTINE collision_interface
|
||||
|
||||
SUBROUTINE findCell_interface(self, part)
|
||||
SUBROUTINE findCell_interface(self, part, oldCell)
|
||||
USE moduleSpecies
|
||||
|
||||
IMPORT:: meshVol
|
||||
CLASS(meshVol), INTENT(in):: self
|
||||
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
|
||||
CLASS(particle), INTENT(inout):: part
|
||||
|
||||
END SUBROUTINE findCell_interface
|
||||
|
|
|
|||
|
|
@ -404,43 +404,47 @@ MODULE moduleMeshCyl
|
|||
tensorS = outerProduct(part%v, part%v)
|
||||
|
||||
vertex => self%n1%output(part%pt)
|
||||
vertex%den = vertex%den + w_p(1)
|
||||
vertex%mom(:) = vertex%mom(:) + w_p(1)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(1)*tensorS
|
||||
vertex%den = vertex%den + part%weight*w_p(1)
|
||||
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
|
||||
|
||||
vertex => self%n2%output(part%pt)
|
||||
vertex%den = vertex%den + w_p(2)
|
||||
vertex%mom(:) = vertex%mom(:) + w_p(2)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(2)*tensorS
|
||||
vertex%den = vertex%den + part%weight*w_p(2)
|
||||
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
|
||||
|
||||
vertex => self%n3%output(part%pt)
|
||||
vertex%den = vertex%den + w_p(3)
|
||||
vertex%mom(:) = vertex%mom(:) + w_p(3)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(3)*tensorS
|
||||
vertex%den = vertex%den + part%weight*w_p(3)
|
||||
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
|
||||
|
||||
vertex => self%n4%output(part%pt)
|
||||
vertex%den = vertex%den + w_p(4)
|
||||
vertex%mom(:) = vertex%mom(:) + w_p(4)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + w_p(4)*tensorS
|
||||
vertex%den = vertex%den + part%weight*w_p(4)
|
||||
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:)
|
||||
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
|
||||
|
||||
|
||||
END SUBROUTINE scatterQuad
|
||||
|
||||
RECURSIVE SUBROUTINE findCellCylQuad(self, part)
|
||||
RECURSIVE SUBROUTINE findCellCylQuad(self, part, oldCell)
|
||||
USE moduleSpecies
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshVolCylQuad), INTENT(in):: self
|
||||
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
|
||||
CLASS(particle), INTENT(inout):: part
|
||||
REAL(8):: xLog(1:2)
|
||||
REAL(8):: xLogArray(1:4)
|
||||
CLASS(*), POINTER:: nextElement
|
||||
INTEGER:: nextInt
|
||||
|
||||
|
||||
xLog = self%phy2log(part%r(1:2))
|
||||
IF (self%inside(xLog(1), xLog(2))) THEN
|
||||
!Checks if particle is inside of current cell
|
||||
IF (PRESENT(oldCell)) THEN
|
||||
!If oldCell, recalculate particle weight, as particle has entered a new cell.
|
||||
part%weight = part%weight*oldCell%volume/self%volume
|
||||
END IF
|
||||
part%e_p = self%n
|
||||
part%xLog = xLog
|
||||
ELSE
|
||||
|
|
@ -464,7 +468,7 @@ MODULE moduleMeshCyl
|
|||
SELECT TYPE(nextElement)
|
||||
CLASS IS(meshVolCyl)
|
||||
!Particle moved to new cell, repeat find procedure
|
||||
CALL nextElement%findCell(part)
|
||||
CALL nextElement%findCell(part, self)
|
||||
|
||||
CLASS IS (meshEdgeCyl)
|
||||
!Particle encountered an edge, execute boundary
|
||||
|
|
|
|||
|
|
@ -60,9 +60,12 @@ MODULE moduleOutput
|
|||
formatValues%velocity = 0.D0
|
||||
formatValues%pressure = 0.D0
|
||||
formatValues%temperature = 0.D0
|
||||
tempVol = speciesIn%weight/(nodeVol*Vol_ref)
|
||||
tempVol = 1.D0/(nodeVol*Vol_ref)
|
||||
IF (rawValues%den > 0.D0) THEN
|
||||
tempVel = rawValues%mom(:)/rawValues%den
|
||||
IF ((tempVel(1) - 1.D0) .EQ. tempVel(1)) THEN
|
||||
PRINT *, rawValues%mom
|
||||
END IF
|
||||
tensorTemp = (rawValues%tensorS(:,:) - rawValues%den*outerProduct(tempVel,tempVel))
|
||||
formatValues%density = rawValues%den*tempVol
|
||||
formatValues%velocity(:) = tempVel
|
||||
|
|
|
|||
|
|
@ -34,6 +34,7 @@ MODULE moduleSpecies
|
|||
INTEGER:: e_p !Index of element in which the particle is located
|
||||
REAL(8):: xLog(1:2) !Logical coordinates of particle in element e_p.
|
||||
LOGICAL:: n_in !Flag that indicates if a particle is in the domain
|
||||
REAL(8):: weight=0.D0
|
||||
|
||||
END TYPE particle
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue