Merge branch 'development' into feature/see
This commit is contained in:
commit
51f2726c3d
23 changed files with 424 additions and 242 deletions
|
|
@ -27,6 +27,10 @@ PROGRAM fpakc
|
|||
|
||||
!Reads the json configuration file
|
||||
CALL readConfig(inputFile)
|
||||
|
||||
!Create output folder and initial files
|
||||
CALL initOutput(inputFile)
|
||||
|
||||
!Do '0' iteration
|
||||
t = tInitial
|
||||
|
||||
|
|
|
|||
|
|
@ -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 = a + FLOOR((b+1-a)*rnd01)
|
||||
|
||||
END FUNCTION randomIntAB
|
||||
|
||||
|
|
@ -91,10 +91,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
|
||||
|
||||
|
|
|
|||
|
|
@ -96,7 +96,7 @@ MODULE moduleTable
|
|||
f = self%fMax
|
||||
|
||||
ELSE
|
||||
i = MINLOC(x - self%x, 1)
|
||||
i = MINLOC(ABS(x - self%x), 1)
|
||||
deltaX = x - self%x(i)
|
||||
IF (deltaX < 0 ) THEN
|
||||
i = i - 1
|
||||
|
|
|
|||
|
|
@ -84,20 +84,6 @@ MODULE moduleInput
|
|||
CALL readParallel(config)
|
||||
CALL checkStatus(config, "readParallel")
|
||||
|
||||
!If everything is correct, creates the output folder
|
||||
CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
|
||||
!Copies input file to output folder
|
||||
CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
|
||||
!Copies particle mesh
|
||||
IF (mesh%dimen > 0) THEN
|
||||
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
|
||||
IF (doubleMesh) THEN
|
||||
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
|
||||
|
||||
END IF
|
||||
|
||||
END IF
|
||||
|
||||
END SUBROUTINE readConfig
|
||||
|
||||
!Checks the status of the JSON case file and, if failed, exits the execution.
|
||||
|
|
@ -322,7 +308,7 @@ MODULE moduleInput
|
|||
LOGICAL:: found
|
||||
CHARACTER(:), ALLOCATABLE:: object
|
||||
INTEGER:: nInitial
|
||||
INTEGER:: i, j, p, e
|
||||
INTEGER:: i, p, e
|
||||
CHARACTER(LEN=2):: iString
|
||||
CHARACTER(:), ALLOCATABLE:: spName
|
||||
INTEGER:: sp
|
||||
|
|
@ -338,7 +324,8 @@ MODULE moduleInput
|
|||
REAL(8):: densityCen
|
||||
!Mean velocity and temperature at particle position
|
||||
REAL(8):: velocityXi(1:3), temperatureXi
|
||||
INTEGER:: nNewPart = 0.D0
|
||||
INTEGER:: nNewPart = 0
|
||||
REAL(8):: weight = 0.D0
|
||||
CLASS(meshCell), POINTER:: cell
|
||||
TYPE(particle), POINTER:: partNew
|
||||
REAL(8):: vTh
|
||||
|
|
@ -357,6 +344,9 @@ MODULE moduleInput
|
|||
!Reads node values at the nodes
|
||||
filename = path // spFile
|
||||
CALL mesh%readInitial(filename, density, velocity, temperature)
|
||||
!Check if initial number of particles is given
|
||||
CALL config%get(object // '.particlesPerCell', nNewPart, found)
|
||||
|
||||
!For each volume in the node, create corresponding particles
|
||||
DO e = 1, mesh%numCells
|
||||
!Scale variables
|
||||
|
|
@ -369,7 +359,11 @@ MODULE moduleInput
|
|||
densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar)
|
||||
|
||||
!Calculate number of particles
|
||||
nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
|
||||
IF (.NOT. found) THEN
|
||||
nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
|
||||
|
||||
END IF
|
||||
weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart)
|
||||
|
||||
!Allocate new particles
|
||||
DO p = 1, nNewPart
|
||||
|
|
@ -406,7 +400,7 @@ MODULE moduleInput
|
|||
|
||||
partNew%n_in = .TRUE.
|
||||
|
||||
partNew%weight = species(sp)%obj%weight
|
||||
partNew%weight = weight
|
||||
|
||||
!Assign particle to temporal list of particles
|
||||
CALL partInitial%add(partNew)
|
||||
|
|
@ -634,7 +628,7 @@ MODULE moduleInput
|
|||
INTEGER:: i, k, ij
|
||||
INTEGER:: pt_i, pt_j
|
||||
REAL(8):: energyThreshold, energyBinding
|
||||
CHARACTER(:), ALLOCATABLE:: electron
|
||||
CHARACTER(:), ALLOCATABLE:: electron, electronSecondary
|
||||
INTEGER:: e
|
||||
CLASS(meshCell), POINTER:: cell
|
||||
|
||||
|
|
@ -711,8 +705,16 @@ MODULE moduleInput
|
|||
IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
|
||||
CALL config%get(object // '.electron', electron, found)
|
||||
IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
|
||||
CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
|
||||
crossSecFilePath, energyThreshold, electron)
|
||||
CALL config%get(object // '.electronSecondary', electronSecondary, found)
|
||||
IF (found) THEN
|
||||
CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
|
||||
crossSecFilePath, energyThreshold, electron, electronSecondary)
|
||||
|
||||
ELSE
|
||||
CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
|
||||
crossSecFilePath, energyThreshold, electron)
|
||||
|
||||
END IF
|
||||
|
||||
CASE ('recombination')
|
||||
!Electorn impact ionization
|
||||
|
|
@ -805,8 +807,8 @@ MODULE moduleInput
|
|||
REAL(8), DIMENSION(:), ALLOCATABLE:: v0
|
||||
REAL(8):: effTime
|
||||
REAL(8):: eThreshold !Energy threshold
|
||||
INTEGER:: speciesID
|
||||
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield
|
||||
INTEGER:: speciesID, electronSecondaryID
|
||||
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield, electronSecondary
|
||||
LOGICAL:: found
|
||||
INTEGER:: nTypes
|
||||
|
||||
|
|
@ -861,8 +863,17 @@ MODULE moduleInput
|
|||
CALL config%get(object // '.crossSection', crossSection, found)
|
||||
IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')
|
||||
|
||||
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
|
||||
speciesID, effTime, crossSection, eThreshold)
|
||||
CALL config%get(object // '.electronSecondary', electronSecondary, found)
|
||||
electronSecondaryID = speciesName2Index(electronSecondary)
|
||||
IF (found) THEN
|
||||
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
|
||||
speciesID, effTime, crossSection, eThreshold,electronSecondaryID)
|
||||
|
||||
ELSE
|
||||
CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
|
||||
speciesID, effTime, crossSection, eThreshold)
|
||||
|
||||
END IF
|
||||
|
||||
CASE('wallTemperature')
|
||||
CALL config%get(object // '.temperature', Tw, found)
|
||||
|
|
@ -921,7 +932,6 @@ MODULE moduleInput
|
|||
LOGICAL:: found
|
||||
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
|
||||
REAL(8):: volume
|
||||
CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT
|
||||
|
||||
object = 'geometry'
|
||||
|
||||
|
|
@ -1234,6 +1244,7 @@ MODULE moduleInput
|
|||
REAL(8):: flow
|
||||
CHARACTER(:), ALLOCATABLE:: units
|
||||
INTEGER:: physicalSurface
|
||||
INTEGER:: particlesPerEdge
|
||||
INTEGER:: sp
|
||||
|
||||
CALL config%info('inject', found, n_children = nInject)
|
||||
|
|
@ -1258,8 +1269,10 @@ MODULE moduleInput
|
|||
CALL config%get(object // '.flow', flow, found)
|
||||
CALL config%get(object // '.units', units, found)
|
||||
CALL config%get(object // '.physicalSurface', physicalSurface, found)
|
||||
particlesPerEdge = 0
|
||||
CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found)
|
||||
|
||||
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
|
||||
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge)
|
||||
|
||||
CALL readVelDistr(config, inject(i), object)
|
||||
|
||||
|
|
@ -1381,5 +1394,37 @@ MODULE moduleInput
|
|||
|
||||
END SUBROUTINE readParallel
|
||||
|
||||
SUBROUTINE initOutput(inputFile)
|
||||
USE moduleRefParam
|
||||
USE moduleMesh, ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl
|
||||
USE moduleOutput, ONLY: path, folder
|
||||
IMPLICIT NONE
|
||||
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
|
||||
INTEGER:: fileReference = 30
|
||||
!If everything is correct, creates the output folder
|
||||
CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
|
||||
!Copies input file to output folder
|
||||
CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
|
||||
!Copies particle mesh
|
||||
IF (mesh%dimen > 0) THEN
|
||||
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
|
||||
IF (doubleMesh) THEN
|
||||
CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)
|
||||
|
||||
END IF
|
||||
|
||||
END IF
|
||||
|
||||
! Write commit of fpakc
|
||||
CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpakc_commit.txt')
|
||||
|
||||
! Write file with reference values
|
||||
OPEN (fileReference, file=path // folder // '/' // 'reference.txt')
|
||||
WRITE(fileReference, "(7(1X,A20))") 'L_ref', 'v_ref', 'ti_ref', 'Vol_ref', 'EF_ref', 'Volt_ref', 'B_ref'
|
||||
WRITE(fileReference, "(7(1X,ES20.6E3))") L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
|
||||
CLOSE(fileReference)
|
||||
|
||||
END SUBROUTINE initOutput
|
||||
|
||||
END MODULE moduleInput
|
||||
|
|
|
|||
|
|
@ -104,6 +104,7 @@ MODULE moduleMesh1DCart
|
|||
USE moduleSpecies
|
||||
USE moduleBoundary
|
||||
USE moduleErrors
|
||||
USE moduleRefParam, ONLY: L_ref
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshEdge1DCart), INTENT(out):: self
|
||||
|
|
@ -122,6 +123,8 @@ MODULE moduleMesh1DCart
|
|||
|
||||
self%x = r1(1)
|
||||
|
||||
self%surface = 1.D0 / L_ref**2
|
||||
|
||||
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
|
||||
|
||||
!Boundary index
|
||||
|
|
|
|||
|
|
@ -104,6 +104,7 @@ MODULE moduleMesh1DRad
|
|||
USE moduleSpecies
|
||||
USE moduleBoundary
|
||||
USE moduleErrors
|
||||
USE moduleRefParam, ONLY: L_ref
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshEdge1DRad), INTENT(out):: self
|
||||
|
|
@ -122,6 +123,8 @@ MODULE moduleMesh1DRad
|
|||
|
||||
self%r = r1(1)
|
||||
|
||||
self%surface = 1.D0 / L_ref**2
|
||||
|
||||
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
|
||||
|
||||
!Boundary index
|
||||
|
|
|
|||
|
|
@ -144,6 +144,7 @@ MODULE moduleMesh2DCart
|
|||
USE moduleSpecies
|
||||
USE moduleBoundary
|
||||
USE moduleErrors
|
||||
USE moduleRefParam, ONLY: L_ref
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshEdge2DCart), INTENT(out):: self
|
||||
|
|
@ -163,7 +164,7 @@ MODULE moduleMesh2DCart
|
|||
r2 = self%n2%getCoordinates()
|
||||
self%x = (/r1(1), r2(1)/)
|
||||
self%y = (/r1(2), r2(2)/)
|
||||
self%weight = 1.D0
|
||||
self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref
|
||||
!Normal vector
|
||||
self%normal = (/ -(self%y(2)-self%y(1)), &
|
||||
self%x(2)-self%x(1) , &
|
||||
|
|
@ -318,6 +319,8 @@ MODULE moduleMesh2DCart
|
|||
INTEGER, INTENT(in):: nNodes
|
||||
REAL(8):: fPsi(1:nNodes)
|
||||
|
||||
fPsi = 0.D0
|
||||
|
||||
fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
|
||||
(1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
|
||||
(1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
|
||||
|
|
@ -508,15 +511,15 @@ MODULE moduleMesh2DCart
|
|||
conv = 1.D0
|
||||
XiO = 0.D0
|
||||
|
||||
f(3) = 0.D0
|
||||
DO WHILE(conv > 1.D-4)
|
||||
dPsi = self%dPsi(XiO, 4)
|
||||
pDer = self%partialDer(4, dPsi)
|
||||
detJ = self%detJac(pDer)
|
||||
invJ = self%invJac(pDer)
|
||||
fPsi = self%fPsi(XiO, 4)
|
||||
f = (/ DOT_PRODUCT(fPsi,self%x), &
|
||||
DOT_PRODUCT(fPsi,self%y), &
|
||||
0.D0 /) - r
|
||||
f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), &
|
||||
DOT_PRODUCT(fPsi,self%y) /) - r(1:2)
|
||||
Xi = XiO - MATMUL(invJ, f)/detJ
|
||||
conv = MAXVAL(DABS(Xi-XiO),1)
|
||||
XiO = Xi
|
||||
|
|
@ -554,6 +557,7 @@ MODULE moduleMesh2DCart
|
|||
|
||||
!Compute element volume
|
||||
PURE SUBROUTINE volumeQuad(self)
|
||||
USE moduleRefParam, ONLY: L_ref
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshCell2DCartQuad), INTENT(inout):: self
|
||||
|
|
@ -569,8 +573,9 @@ MODULE moduleMesh2DCart
|
|||
pDer = self%partialDer(4, dPsi)
|
||||
detJ = self%detJac(pDer)
|
||||
fPsi = self%fPsi(Xi, 4)
|
||||
|
||||
!Compute total volume of the cell
|
||||
self%volume = detJ*4.D0
|
||||
self%volume = detJ*4.D0/L_ref
|
||||
!Compute volume per node
|
||||
self%n1%v = self%n1%v + fPsi(1)*self%volume
|
||||
self%n2%v = self%n2%v + fPsi(2)*self%volume
|
||||
|
|
@ -762,6 +767,7 @@ MODULE moduleMesh2DCart
|
|||
pDer = self%partialDer(3, dPsi)
|
||||
detJ = self%detJac(pDer)
|
||||
invJ = self%invJac(pDer)
|
||||
|
||||
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
|
||||
|
||||
END DO
|
||||
|
|
|
|||
|
|
@ -144,6 +144,7 @@ MODULE moduleMesh2DCyl
|
|||
USE moduleSpecies
|
||||
USE moduleBoundary
|
||||
USE moduleErrors
|
||||
USE moduleConstParam, ONLY: PI
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshEdge2DCyl), INTENT(out):: self
|
||||
|
|
@ -163,7 +164,15 @@ 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
|
||||
!Edge surface
|
||||
IF (self%z(2) /= self%z(1)) THEN
|
||||
self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1))
|
||||
|
||||
ELSE
|
||||
self%surface = ABS(self%r(2)**2 - self%r(1)**2)
|
||||
|
||||
END IF
|
||||
self%surface = self%surface * PI
|
||||
!Normal vector
|
||||
self%normal = (/ -(self%r(2)-self%r(1)), &
|
||||
self%z(2)-self%z(1) , &
|
||||
|
|
@ -223,21 +232,13 @@ 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
|
||||
|
|
@ -246,7 +247,6 @@ MODULE moduleMesh2DCyl
|
|||
!QUAD FUNCTIONS
|
||||
!Init element
|
||||
SUBROUTINE initCellQuad2DCyl(self, n, p, nodes)
|
||||
USE moduleRefParam
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshCell2DCylQuad), INTENT(out):: self
|
||||
|
|
@ -326,6 +326,8 @@ MODULE moduleMesh2DCyl
|
|||
INTEGER, INTENT(in):: nNodes
|
||||
REAL(8):: fPsi(1:nNodes)
|
||||
|
||||
fPsi = 0.D0
|
||||
|
||||
fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
|
||||
(1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
|
||||
(1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
|
||||
|
|
@ -496,7 +498,7 @@ MODULE moduleMesh2DCyl
|
|||
|
||||
END FUNCTION elemFQuad
|
||||
|
||||
!Checks if Xi is inside the element
|
||||
!Check if Xi is inside the element
|
||||
PURE FUNCTION insideQuad(Xi) RESULT(ins)
|
||||
IMPLICIT NONE
|
||||
|
||||
|
|
@ -524,15 +526,15 @@ MODULE moduleMesh2DCyl
|
|||
conv = 1.D0
|
||||
XiO = 0.D0
|
||||
|
||||
f(3) = 0.D0
|
||||
DO WHILE(conv > 1.D-4)
|
||||
dPsi = self%dPsi(XiO, 4)
|
||||
pDer = self%partialDer(4, dPsi)
|
||||
detJ = self%detJac(pDer)
|
||||
invJ = self%invJac(pDer)
|
||||
fPsi = self%fPsi(XiO, 4)
|
||||
f = (/ DOT_PRODUCT(fPsi,self%z), &
|
||||
DOT_PRODUCT(fPsi,self%r), &
|
||||
0.D0 /) - r
|
||||
f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), &
|
||||
DOT_PRODUCT(fPsi,self%r) /) - r(1:2)
|
||||
Xi = XiO - MATMUL(invJ, f)/detJ
|
||||
conv = MAXVAL(DABS(Xi-XiO),1)
|
||||
XiO = Xi
|
||||
|
|
@ -553,7 +555,7 @@ MODULE moduleMesh2DCyl
|
|||
|
||||
XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
|
||||
nextInt = MAXLOC(XiArray,1)
|
||||
!Selects the higher value of directions and searches in that direction
|
||||
!Select the higher value of directions and searches in that direction
|
||||
NULLIFY(neighbourElement)
|
||||
SELECT CASE (nextInt)
|
||||
CASE (1)
|
||||
|
|
@ -581,6 +583,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)
|
||||
|
|
@ -589,18 +592,18 @@ MODULE moduleMesh2DCyl
|
|||
fPsi = self%fPsi(Xi, 4)
|
||||
r = DOT_PRODUCT(fPsi,self%r)
|
||||
!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/)
|
||||
self%volume = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral)
|
||||
!Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis.
|
||||
Xi = (/-5.D-1, -0.25D0, 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/)
|
||||
Xi = (/ 5.D-1, -0.25D0, 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/)
|
||||
Xi = (/ 5.D-1, 0.75D0, 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/)
|
||||
Xi = (/-5.D-1, 0.75D0, 0.D0/)
|
||||
r = self%gatherF(Xi, 4, self%r)
|
||||
self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8
|
||||
|
||||
|
|
|
|||
|
|
@ -109,6 +109,7 @@ MODULE moduleMesh3DCart
|
|||
USE moduleBoundary
|
||||
USE moduleErrors
|
||||
USE moduleMath
|
||||
USE moduleRefParam, ONLY: L_ref
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshEdge3DCartTria), INTENT(out):: self
|
||||
|
|
@ -142,6 +143,8 @@ MODULE moduleMesh3DCart
|
|||
self%normal = crossProduct(vec1, vec2)
|
||||
self%normal = normalize(self%normal)
|
||||
|
||||
self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D
|
||||
|
||||
!Boundary index
|
||||
self%boundary => boundary(bt)
|
||||
ALLOCATE(self%fBoundary(1:nSpecies))
|
||||
|
|
|
|||
|
|
@ -108,6 +108,7 @@ MODULE moduleMeshInputGmsh2
|
|||
READ(10, *) totalNumElem
|
||||
|
||||
!Count edges and volume elements
|
||||
numEdges = 0
|
||||
SELECT TYPE(self)
|
||||
TYPE IS(meshParticles)
|
||||
self%numEdges = 0
|
||||
|
|
@ -328,7 +329,7 @@ MODULE moduleMeshInputGmsh2
|
|||
|
||||
DO i = 1, numNodes
|
||||
!Reads the density
|
||||
READ(10, *), e, density(i)
|
||||
READ(10, *) e, density(i)
|
||||
|
||||
END DO
|
||||
|
||||
|
|
@ -339,7 +340,7 @@ MODULE moduleMeshInputGmsh2
|
|||
|
||||
DO i = 1, numNodes
|
||||
!Reads the velocity
|
||||
READ(10, *), e, velocity(i, 1:3)
|
||||
READ(10, *) e, velocity(i, 1:3)
|
||||
|
||||
END DO
|
||||
|
||||
|
|
|
|||
|
|
@ -275,6 +275,7 @@ MODULE moduleMeshInputVTU
|
|||
END DO
|
||||
|
||||
!Count the number of edges
|
||||
numEdges = 0
|
||||
SELECT CASE(self%dimen)
|
||||
CASE(3)
|
||||
!Edges are triangles, type 5 in VTK
|
||||
|
|
@ -495,7 +496,7 @@ MODULE moduleMeshInputVTU
|
|||
END SELECT
|
||||
|
||||
END DO
|
||||
|
||||
|
||||
!Call mesh connectivity
|
||||
CALL self%connectMesh
|
||||
|
||||
|
|
|
|||
|
|
@ -209,7 +209,7 @@ MODULE moduleMeshOutputVTU
|
|||
WRITE(fileID,"(8X,A)") '<CellData>'
|
||||
!Electric field
|
||||
WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
|
||||
WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
|
||||
WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
|
||||
WRITE(fileID,"(10X,A)") '</DataArray>'
|
||||
WRITE(fileID,"(8X,A)") '</CellData>'
|
||||
|
||||
|
|
@ -315,9 +315,8 @@ MODULE moduleMeshOutputVTU
|
|||
|
||||
CLASS(meshParticles), INTENT(in):: self
|
||||
INTEGER, INTENT(in):: t
|
||||
INTEGER:: n, i, fileID
|
||||
INTEGER:: i, fileID
|
||||
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
|
||||
TYPE(outputFormat):: output(1:self%numNodes)
|
||||
|
||||
fileID = 60
|
||||
|
||||
|
|
@ -352,10 +351,9 @@ MODULE moduleMeshOutputVTU
|
|||
|
||||
CLASS(meshGeneric), INTENT(in):: self
|
||||
INTEGER, INTENT(in):: t
|
||||
INTEGER:: n, i, fileID
|
||||
INTEGER:: fileID
|
||||
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
|
||||
CHARACTER (LEN=iterationDigits):: tstring
|
||||
TYPE(outputFormat):: output(1:self%numNodes)
|
||||
|
||||
fileID = 62
|
||||
|
||||
|
|
@ -424,9 +422,8 @@ MODULE moduleMeshOutputVTU
|
|||
IMPLICIT NONE
|
||||
|
||||
CLASS(meshParticles), INTENT(in):: self
|
||||
INTEGER:: n, i, fileIDMean, fileIDDeviation
|
||||
INTEGER:: i, fileIDMean, fileIDDeviation
|
||||
CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation
|
||||
TYPE(outputFormat):: output(1:self%numNodes)
|
||||
|
||||
fileIDMean = 66
|
||||
fileIDDeviation = 67
|
||||
|
|
|
|||
|
|
@ -76,8 +76,8 @@ MODULE moduleMesh
|
|||
CLASS(meshCell), POINTER:: eColl => NULL()
|
||||
!Normal vector
|
||||
REAL(8):: normal(1:3)
|
||||
!Weight for random injection of particles
|
||||
REAL(8):: weight = 1.D0
|
||||
! Surface of edge
|
||||
REAL(8):: surface = 0.D0
|
||||
!Pointer to boundary type
|
||||
TYPE(boundaryCont), POINTER:: boundary
|
||||
!Array of functions for boundary conditions
|
||||
|
|
@ -613,6 +613,7 @@ MODULE moduleMesh
|
|||
INTEGER:: sp
|
||||
INTEGER:: i
|
||||
CLASS(meshNode), POINTER:: node
|
||||
REAL(8):: pFraction !Particle fraction
|
||||
|
||||
cellNodes = self%getNodes(nNodes)
|
||||
fPsi = self%fPsi(part%Xi, nNodes)
|
||||
|
|
@ -623,10 +624,11 @@ MODULE moduleMesh
|
|||
|
||||
DO i = 1, nNodes
|
||||
node => mesh%nodes(cellNodes(i))%obj
|
||||
pFraction = fPsi(i)*part%weight
|
||||
CALL OMP_SET_LOCK(node%lock)
|
||||
node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i)
|
||||
node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:)
|
||||
node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS
|
||||
node%output(sp)%den = node%output(sp)%den + pFraction
|
||||
node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:)
|
||||
node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
|
||||
CALL OMP_UNSET_LOCK(node%lock)
|
||||
|
||||
END DO
|
||||
|
|
@ -911,7 +913,9 @@ MODULE moduleMesh
|
|||
!Loop over collisions
|
||||
DO c = 1, interactionMatrix(k)%amount
|
||||
IF (rnd_real <= probabilityColl(c)) THEN
|
||||
!$OMP CRITICAL
|
||||
CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
|
||||
!$OMP END CRITICAL
|
||||
|
||||
!If collisions are gonna be output, count the collision
|
||||
IF (collOutput) THEN
|
||||
|
|
@ -1021,6 +1025,9 @@ MODULE moduleMesh
|
|||
ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3))
|
||||
ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3))
|
||||
ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount))
|
||||
deltaV_ij = 0.D0
|
||||
p_ij = 0.D0
|
||||
mass_ij = 0.D0
|
||||
!Loop over particles of species_i
|
||||
partTemp => cell%listPart_in(i)%head
|
||||
p = 1
|
||||
|
|
@ -1105,6 +1112,9 @@ MODULE moduleMesh
|
|||
ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3))
|
||||
ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3))
|
||||
ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount))
|
||||
deltaV_ji = 0.D0
|
||||
p_ji = 0.D0
|
||||
mass_ji = 0.D0
|
||||
!Loop over particles of species_j
|
||||
partTemp => cell%listPart_in(j)%head
|
||||
p = 1
|
||||
|
|
|
|||
|
|
@ -147,7 +147,13 @@ MODULE moduleMeshBoundary
|
|||
ALLOCATE(newElectron)
|
||||
ALLOCATE(newIon)
|
||||
|
||||
newElectron%species => part%species
|
||||
IF (ASSOCIATED(bound%electronSecondary)) THEN
|
||||
newElectron%species => bound%electronSecondary
|
||||
|
||||
ELSE
|
||||
newElectron%species => part%species
|
||||
|
||||
END IF
|
||||
newIon%species => bound%species
|
||||
|
||||
newElectron%v = v0 + (1.D0 + bound%deltaV*v0/NORM2(v0))
|
||||
|
|
|
|||
|
|
@ -38,6 +38,7 @@ MODULE moduleBoundary
|
|||
TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
|
||||
REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
|
||||
CLASS(speciesGeneric), POINTER:: species !Ion species
|
||||
CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
|
||||
TYPE(table1D):: crossSection
|
||||
REAL(8):: effectiveTime
|
||||
REAL(8):: eThreshold
|
||||
|
|
@ -113,17 +114,19 @@ MODULE moduleBoundary
|
|||
|
||||
END SUBROUTINE initWallTemperature
|
||||
|
||||
SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, speciesID, effTime, crossSection, eThreshold)
|
||||
SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
|
||||
USE moduleRefParam
|
||||
USE moduleSpecies
|
||||
USE moduleCaseParam
|
||||
USE moduleConstParam
|
||||
USE moduleErrors
|
||||
IMPLICIT NONE
|
||||
|
||||
CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
|
||||
REAL(8), INTENT(in):: me !Electron mass
|
||||
REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
|
||||
INTEGER:: speciesID
|
||||
INTEGER, INTENT(in):: ion
|
||||
INTEGER, OPTIONAL, INTENT(in):: electronSecondary
|
||||
REAL(8):: effTime
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
|
||||
REAL(8), INTENT(in):: eThreshold
|
||||
|
|
@ -136,7 +139,22 @@ MODULE moduleBoundary
|
|||
boundary%n0 = n0 * Vol_ref
|
||||
boundary%v0 = v0 / v_ref
|
||||
boundary%vTh = DSQRT(kb*T0/m0)/v_ref
|
||||
boundary%species => species(speciesID)%obj
|
||||
boundary%species => species(ion)%obj
|
||||
IF (PRESENT(electronSecondary)) THEN
|
||||
SELECT TYPE(sp => species(electronSecondary)%obj)
|
||||
TYPE IS(speciesCharged)
|
||||
boundary%electronSecondary => sp
|
||||
|
||||
CLASS DEFAULT
|
||||
CALL criticalError("Species " // sp%name // " chosen for " // &
|
||||
"secondary electron is not a charged species", 'initIonization')
|
||||
|
||||
END SELECT
|
||||
|
||||
ELSE
|
||||
boundary%electronSecondary => NULL()
|
||||
|
||||
END IF
|
||||
boundary%effectiveTime = effTime / ti_ref
|
||||
CALL boundary%crossSection%init(crossSection)
|
||||
CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
|
||||
|
|
|
|||
|
|
@ -43,7 +43,8 @@ MODULE moduleCollisions
|
|||
TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization
|
||||
REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization
|
||||
REAL(8):: deltaV !Change in velocity due to exchange of eThreshold
|
||||
CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
|
||||
CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
|
||||
CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: collide => collideBinaryIonization
|
||||
|
||||
|
|
@ -241,7 +242,7 @@ MODULE moduleCollisions
|
|||
|
||||
!ELECTRON IMPACT IONIZATION
|
||||
!Inits electron impact ionization
|
||||
SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron)
|
||||
SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, electron, electronSecondary)
|
||||
USE moduleTable
|
||||
USE moduleRefParam
|
||||
USE moduleConstParam
|
||||
|
|
@ -253,7 +254,8 @@ MODULE moduleCollisions
|
|||
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
|
||||
REAL(8), INTENT(in):: energyThreshold
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
|
||||
INTEGER:: electronIndex
|
||||
CHARACTER(:), ALLOCATABLE, OPTIONAL, INTENT(in):: electronSecondary
|
||||
INTEGER:: electronIndex, electronSecondaryIndex
|
||||
|
||||
ALLOCATE(collisionBinaryIonization:: collision)
|
||||
|
||||
|
|
@ -278,10 +280,27 @@ MODULE moduleCollisions
|
|||
|
||||
CLASS DEFAULT
|
||||
CALL criticalError("Species " // sp%name // " chosen for " // &
|
||||
"secondary electron is not a charged species", 'initBinaryIonization')
|
||||
"impacting electron is not a charged species", 'initBinaryIonization')
|
||||
|
||||
END SELECT
|
||||
|
||||
IF (PRESENT(electronSecondary)) THEN
|
||||
electronSecondaryIndex = speciesName2Index(electronSecondary)
|
||||
SELECT TYPE(sp => species(electronSecondaryIndex)%obj)
|
||||
TYPE IS(speciesCharged)
|
||||
collision%electronSecondary => sp
|
||||
|
||||
CLASS DEFAULT
|
||||
CALL criticalError("Species " // sp%name // " chosen for " // &
|
||||
"secondary electron is not a charged species", 'initBinaryIonization')
|
||||
|
||||
END SELECT
|
||||
|
||||
ELSE
|
||||
collision%electronSecondary => NULL()
|
||||
|
||||
END IF
|
||||
|
||||
!momentum change per ionization process
|
||||
collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
|
||||
|
||||
|
|
@ -336,6 +355,12 @@ MODULE moduleCollisions
|
|||
!Copy basic information from primary electron
|
||||
newElectron = electron
|
||||
|
||||
!If secondary electron species indicates, convert
|
||||
IF (ASSOCIATED(self%electronSecondary)) THEN
|
||||
newElectron%species => self%electronSecondary
|
||||
|
||||
END IF
|
||||
|
||||
!Secondary electorn gains energy from ionization
|
||||
newElectron%v = vChange
|
||||
|
||||
|
|
@ -362,7 +387,7 @@ MODULE moduleCollisions
|
|||
CALL sp%ionize(neutral)
|
||||
|
||||
CLASS DEFAULT
|
||||
! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
|
||||
CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
|
||||
RETURN
|
||||
|
||||
END SELECT
|
||||
|
|
|
|||
|
|
@ -61,8 +61,9 @@ MODULE moduleInject
|
|||
CLASS(speciesGeneric), POINTER:: species !Species of injection
|
||||
INTEGER:: nEdges
|
||||
INTEGER, ALLOCATABLE:: edges(:) !Array with edges
|
||||
REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
|
||||
REAL(8):: sumWeight
|
||||
INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge
|
||||
REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
|
||||
REAL(8):: surface ! Total surface of injection
|
||||
TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
|
||||
CONTAINS
|
||||
PROCEDURE, PASS:: init => initInject
|
||||
|
|
@ -75,7 +76,7 @@ MODULE moduleInject
|
|||
|
||||
CONTAINS
|
||||
!Initialize an injection of particles
|
||||
SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
|
||||
SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge)
|
||||
USE moduleMesh
|
||||
USE moduleRefParam
|
||||
USE moduleConstParam
|
||||
|
|
@ -87,48 +88,28 @@ MODULE moduleInject
|
|||
CLASS(injectGeneric), INTENT(inout):: self
|
||||
INTEGER, INTENT(in):: i
|
||||
REAL(8), INTENT(in):: v, n(1:3), T(1:3)
|
||||
INTEGER, INTENT(in):: sp, physicalSurface
|
||||
INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
|
||||
REAL(8):: tauInject
|
||||
REAL(8), INTENT(in):: flow
|
||||
CHARACTER(:), ALLOCATABLE, INTENT(in):: units
|
||||
INTEGER:: e, et
|
||||
INTEGER:: phSurface(1:mesh%numEdges)
|
||||
INTEGER:: nVolColl
|
||||
REAL(8):: fluxPerStep = 0.D0
|
||||
|
||||
self%id = i
|
||||
self%vMod = v / v_ref
|
||||
self%n = n / NORM2(n)
|
||||
self%T = T / T_ref
|
||||
self%species => species(sp)%obj
|
||||
tauInject = tau(self%species%n)
|
||||
SELECT CASE(units)
|
||||
CASE ("sccm")
|
||||
!Standard cubic centimeter per minute
|
||||
self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight)
|
||||
|
||||
CASE ("A")
|
||||
!Input current in Ampers
|
||||
self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
|
||||
|
||||
CASE ("part/s")
|
||||
!Input current in Ampers
|
||||
self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
|
||||
|
||||
CASE DEFAULT
|
||||
CALL criticalError("No support for units: " // units, 'initInject')
|
||||
|
||||
END SELECT
|
||||
!Scale particles for different species steps
|
||||
IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
|
||||
|
||||
self%id = i
|
||||
self%vMod = v / v_ref
|
||||
self%n = n / NORM2(n)
|
||||
self%T = T / T_ref
|
||||
!Gets the edge elements from which particles are injected
|
||||
DO e = 1, mesh%numEdges
|
||||
phSurface(e) = mesh%edges(e)%obj%physicalSurface
|
||||
|
||||
END DO
|
||||
|
||||
self%nEdges = COUNT(phSurface == physicalSurface)
|
||||
ALLOCATE(inject(i)%edges(1:self%nEdges))
|
||||
ALLOCATE(self%edges(1:self%nEdges))
|
||||
ALLOCATE(self%particlesPerEdge(1:self%nEdges))
|
||||
ALLOCATE(self%weightPerEdge(1:self%nEdges))
|
||||
et = 0
|
||||
DO e=1, mesh%numEdges
|
||||
IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
|
||||
|
|
@ -160,15 +141,63 @@ MODULE moduleInject
|
|||
|
||||
END DO
|
||||
|
||||
!Calculates cumulative probability
|
||||
ALLOCATE(self%cumWeight(1:self%nEdges))
|
||||
et = 1
|
||||
self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
|
||||
DO et = 2, self%nEdges
|
||||
self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
|
||||
!Calculates total area
|
||||
self%surface = 0.D0
|
||||
DO et = 1, self%nEdges
|
||||
self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface
|
||||
|
||||
END DO
|
||||
self%sumWeight = self%cumWeight(self%nEdges)
|
||||
|
||||
! Information about species and flux
|
||||
self%species => species(sp)%obj
|
||||
tauInject = tau(self%species%n)
|
||||
! Convert units
|
||||
SELECT CASE(units)
|
||||
CASE ("sccm")
|
||||
!Standard cubic centimeter per minute
|
||||
fluxPerStep = flow*sccm2atomPerS
|
||||
|
||||
CASE ("A")
|
||||
!Current in Ampers
|
||||
fluxPerStep = flow/qe
|
||||
|
||||
CASE ("Am2")
|
||||
!Input current in Ampers per square meter
|
||||
fluxPerStep = flow*self%surface*L_ref**2/qe
|
||||
|
||||
CASE ("part/s")
|
||||
!Input current in Ampers
|
||||
fluxPerStep = flow
|
||||
|
||||
CASE DEFAULT
|
||||
CALL criticalError("No support for units: " // units, 'initInject')
|
||||
|
||||
END SELECT
|
||||
fluxPerStep = fluxPerStep * tauInject * ti_ref / self%surface
|
||||
|
||||
!Assign particles per edge
|
||||
IF (particlesPerEdge > 0) THEN
|
||||
! Particles per edge defined by the user
|
||||
self%particlesPerEdge = particlesPerEdge
|
||||
DO et = 1, self%nEdges
|
||||
self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge)
|
||||
|
||||
END DO
|
||||
|
||||
ELSE
|
||||
! No particles assigned per edge, use the species weight
|
||||
self%weightPerEdge = self%species%weight
|
||||
DO et = 1, self%nEdges
|
||||
self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight)
|
||||
|
||||
END DO
|
||||
|
||||
END IF
|
||||
|
||||
self%nParticles = SUM(self%particlesPerEdge)
|
||||
|
||||
!Scale particles for different species steps
|
||||
IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
|
||||
|
||||
END SUBROUTINE initInject
|
||||
|
||||
|
|
@ -279,9 +308,8 @@ MODULE moduleInject
|
|||
IMPLICIT NONE
|
||||
|
||||
CLASS(injectGeneric), INTENT(in):: self
|
||||
INTEGER:: randomX
|
||||
INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
|
||||
INTEGER:: i
|
||||
INTEGER, SAVE:: nMin
|
||||
INTEGER:: i, j, e
|
||||
INTEGER:: n, sp
|
||||
CLASS(meshEdge), POINTER:: randomEdge
|
||||
REAL(8):: direction(1:3)
|
||||
|
|
@ -296,61 +324,66 @@ MODULE moduleInject
|
|||
END IF
|
||||
|
||||
END DO
|
||||
|
||||
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
|
||||
|
||||
!$OMP DO
|
||||
DO n = nMin, nMax
|
||||
randomX = randomWeighted(self%cumWeight, self%sumWeight)
|
||||
DO e = 1, self%nEdges
|
||||
! Select edge for injection
|
||||
randomEdge => mesh%edges(self%edges(e))%obj
|
||||
! Inject particles in edge
|
||||
DO i = 1, self%particlesPerEdge(e)
|
||||
! Index in the global partInj array
|
||||
n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i
|
||||
!Particle is considered to be outside the domain
|
||||
partInj(n)%n_in = .FALSE.
|
||||
!Random position in edge
|
||||
partInj(n)%r = randomEdge%randPos()
|
||||
!Assign weight to particle.
|
||||
partInj(n)%weight = self%weightPerEdge(e)
|
||||
!Volume associated to the edge:
|
||||
IF (ASSOCIATED(randomEdge%e1)) THEN
|
||||
partInj(n)%cell = randomEdge%e1%n
|
||||
|
||||
randomEdge => mesh%edges(self%edges(randomX))%obj
|
||||
!Random position in edge
|
||||
partInj(n)%r = randomEdge%randPos()
|
||||
!Volume associated to the edge:
|
||||
IF (ASSOCIATED(randomEdge%e1)) THEN
|
||||
partInj(n)%cell = randomEdge%e1%n
|
||||
ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
|
||||
partInj(n)%cell = randomEdge%e2%n
|
||||
|
||||
ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
|
||||
partInj(n)%cell = randomEdge%e2%n
|
||||
ELSE
|
||||
CALL criticalError("No Volume associated to edge", 'addParticles')
|
||||
|
||||
ELSE
|
||||
CALL criticalError("No Volume associated to edge", 'addParticles')
|
||||
END IF
|
||||
partInj(n)%cellColl = randomEdge%eColl%n
|
||||
sp = self%species%n
|
||||
|
||||
END IF
|
||||
partInj(n)%cellColl = randomEdge%eColl%n
|
||||
sp = self%species%n
|
||||
!Assign particle type
|
||||
partInj(n)%species => self%species
|
||||
|
||||
!Assign particle type
|
||||
partInj(n)%species => self%species
|
||||
direction = self%n
|
||||
|
||||
direction = self%n
|
||||
partInj(n)%v = 0.D0
|
||||
|
||||
!Sample initial velocity
|
||||
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
|
||||
self%v(2)%obj%randomVel(), &
|
||||
self%v(3)%obj%randomVel() /)
|
||||
!Sample initial velocity
|
||||
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
|
||||
self%v(2)%obj%randomVel(), &
|
||||
self%v(3)%obj%randomVel() /)
|
||||
|
||||
!For each direction, velocities have to agree with the direction of injection
|
||||
DO i = 1, 3
|
||||
DO WHILE (partInj(n)%v(i)*direction(i) < 0)
|
||||
partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel()
|
||||
!For each direction, velocities have to agree with the direction of injection
|
||||
DO j = 1, 3
|
||||
DO WHILE (partInj(n)%v(i)*direction(i) < 0)
|
||||
partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel()
|
||||
|
||||
END DO
|
||||
|
||||
END DO
|
||||
|
||||
END DO
|
||||
!Obtain natural coordinates of particle in cell
|
||||
partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
|
||||
!Push new particle with the minimum time step
|
||||
CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
|
||||
!Assign cell to new particle
|
||||
CALL solver%updateParticleCell(partInj(n))
|
||||
|
||||
!Obtain natural coordinates of particle in cell
|
||||
partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
|
||||
!Push new particle with the minimum time step
|
||||
CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
|
||||
!Assign cell to new particle
|
||||
CALL solver%updateParticleCell(partInj(n))
|
||||
END DO
|
||||
|
||||
END DO
|
||||
!$OMP END DO
|
||||
|
|
|
|||
|
|
@ -101,7 +101,7 @@ MODULE moduleProbe
|
|||
|
||||
!Maximum radius
|
||||
!TODO: Make this an input parameter
|
||||
self%maxR = 1.D0
|
||||
self%maxR = 1.D-2/L_ref
|
||||
|
||||
!Init the probe lock
|
||||
CALL OMP_INIT_LOCK(self%lock)
|
||||
|
|
@ -148,7 +148,7 @@ MODULE moduleProbe
|
|||
deltaR = NORM2(self%r - part%r)
|
||||
|
||||
!Only include particle if it is inside the maximum radius
|
||||
IF (deltaR < self%maxR) THEN
|
||||
! IF (deltaR < self%maxR) THEN
|
||||
!find lower index for all dimensions
|
||||
CALL self%findLowerIndex(part%v, i, j, k, inside)
|
||||
|
||||
|
|
@ -162,28 +162,28 @@ MODULE moduleProbe
|
|||
fk = self%vk(k+1) - part%v(3)
|
||||
fk1 = part%v(3) - self%vk(k)
|
||||
|
||||
! weight = part%weight * DEXP(deltaR/self%maxR)
|
||||
weight = part%weight
|
||||
weight = part%weight * DEXP(-deltaR/self%maxR)
|
||||
! weight = part%weight
|
||||
|
||||
!Lock the probe
|
||||
CALL OMP_SET_LOCK(self%lock)
|
||||
|
||||
!Assign particle weight to distribution function
|
||||
self%f(i , j , k ) = fi * fj * fk * weight
|
||||
self%f(i+1, j , k ) = fi1 * fj * fk * weight
|
||||
self%f(i , j+1, k ) = fi * fj1 * fk * weight
|
||||
self%f(i+1, j+1, k ) = fi1 * fj1 * fk * weight
|
||||
self%f(i , j , k+1) = fi * fj * fk1 * weight
|
||||
self%f(i+1, j , k+1) = fi1 * fj * fk1 * weight
|
||||
self%f(i , j+1, k+1) = fi * fj1 * fk1 * weight
|
||||
self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight
|
||||
self%f(i , j , k ) = self%f(i , j , k ) + fi * fj * fk * weight
|
||||
self%f(i+1, j , k ) = self%f(i+1, j , k ) + fi1 * fj * fk * weight
|
||||
self%f(i , j+1, k ) = self%f(i , j+1, k ) + fi * fj1 * fk * weight
|
||||
self%f(i+1, j+1, k ) = self%f(i+1, j+1, k ) + fi1 * fj1 * fk * weight
|
||||
self%f(i , j , k+1) = self%f(i , j , k+1) + fi * fj * fk1 * weight
|
||||
self%f(i+1, j , k+1) = self%f(i+1, j , k+1) + fi1 * fj * fk1 * weight
|
||||
self%f(i , j+1, k+1) = self%f(i , j+1, k+1) + fi * fj1 * fk1 * weight
|
||||
self%f(i+1, j+1, k+1) = self%f(i+1, j+1, k+1) + fi1 * fj1 * fk1 * weight
|
||||
|
||||
!Unlock the probe
|
||||
CALL OMP_UNSET_LOCK(self%lock)
|
||||
|
||||
END IF
|
||||
|
||||
END IF
|
||||
! END IF
|
||||
|
||||
END IF
|
||||
|
||||
|
|
|
|||
|
|
@ -30,8 +30,9 @@ MODULE moduleEM
|
|||
INTEGER, ALLOCATABLE:: nodes(:)
|
||||
INTEGER:: n
|
||||
|
||||
nNodes = 1
|
||||
nNodes = edge%nNodes
|
||||
nodes = edge%getNodes(nNodes)
|
||||
nodes = edge%getNodes(nNodes)
|
||||
|
||||
DO n = 1, nNodes
|
||||
SELECT CASE(self%typeEM)
|
||||
|
|
|
|||
|
|
@ -517,11 +517,6 @@ MODULE moduleSolver
|
|||
|
||||
INTEGER, INTENT(in):: t
|
||||
|
||||
IF (t == tInitial) THEN
|
||||
CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpack_commit.txt')
|
||||
|
||||
END IF
|
||||
|
||||
CALL outputProbes(t)
|
||||
|
||||
counterOutput = counterOutput + 1
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue