JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'development' into feature/BoltzmannElectrons

Browse source
This commit is contained in:
Jorge Gonzalez 2024-10-13 14:54:34 +02:00
commit 221de46734
26 changed files with 794 additions and 431 deletions
Download patch file Download diff file Expand all files Collapse all files

29
src/fpakc.f90
View file

@ -11,12 +11,12 @@ PROGRAM fpakc
USE OMP_LIB
IMPLICIT NONE
! t = time step
INTEGER:: t
! arg1 = Input argument 1 (input file)
CHARACTER(200):: arg1
! inputFile = path+name of input file
CHARACTER(:), ALLOCATABLE:: inputFile
! generic integer for time step
INTEGER:: t
tStep = omp_get_wtime()
!Gets the input file
@ -27,11 +27,18 @@ PROGRAM fpakc
!Reads the json configuration file
CALL readConfig(inputFile)
!Create output folder and initial files
CALL initOutput(inputFile)
!Do '0' iteration
t = tInitial
timeStep = tInitial
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP SINGLE
! Initial reset of probes
CALL resetProbes()
CALL verboseError("Initial scatter of particles...")
!$OMP END SINGLE
CALL doScatter()
@ -45,19 +52,21 @@ PROGRAM fpakc
tStep = omp_get_wtime() - tStep
!Output initial state
CALL doOutput(t)
CALL doOutput()
CALL verboseError('Starting main loop...')
!$OMP PARALLEL DEFAULT(SHARED)
DO t = tInitial + 1, tFinal
!Insert new particles and push them
!$OMP SINGLE
tStep = omp_get_wtime()
! Update global time step index
timeStep = t
!Checks if a species needs to me moved in this iteration
CALL solver%updatePushSpecies(t)
CALL solver%updatePushSpecies()
!Checks if probes need to be calculated this iteration
CALL resetProbes(t)
CALL resetProbes()
tPush = omp_get_wtime()
!$OMP END SINGLE
@ -75,7 +84,7 @@ PROGRAM fpakc
!$OMP END SINGLE
IF (doMCCollisions) THEN
CALL meshForMCC%doCollisions(t)
CALL meshForMCC%doCollisions()
END IF
@ -120,12 +129,12 @@ PROGRAM fpakc
!$OMP SINGLE
tEMField = omp_get_wtime() - tEMField
CALL doAverage(t)
CALL doAverage()
tStep = omp_get_wtime() - tStep
!Output data
CALL doOutput(t)
CALL doOutput()
!$OMP END SINGLE
END DO

5
src/modules/common/moduleCaseParam.f90
View file

@ -2,8 +2,13 @@
MODULE moduleCaseParam
!Final and initial iterations
INTEGER:: tFinal, tInitial = 0
! Global index of current iteration
INTEGER:: timeStep
! Time step for all species
REAL(8), ALLOCATABLE:: tau(:)
! Minimum time step
REAL(8):: tauMin
! Time step for Monte-Carlo Collisions
REAL(8):: tauColl
END MODULE moduleCaseParam

21
src/modules/common/moduleRandom.f90
View file

@ -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
@ -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

236
src/modules/init/moduleInput.f90
View file

@ -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.
@ -282,8 +268,8 @@ MODULE moduleInput
CALL readEMBoundary(config)
!Read constant magnetic field
DO i = 1, 3
WRITE(istring, '(i2)') i
CALL config%get(object // '.B(' // istring // ')', B(i), found)
WRITE(iString, '(i2)') i
CALL config%get(object // '.B(' // iString // ')', B(i), found)
IF (.NOT. found) THEN
CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')
@ -326,7 +312,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
@ -342,7 +328,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
@ -361,6 +348,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
@ -373,7 +363,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
@ -410,7 +404,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)
@ -809,7 +803,7 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
INTEGER:: i, s
CHARACTER(2):: istring, sString
CHARACTER(2):: iString, sString
CHARACTER(:), ALLOCATABLE:: object, bType
REAL(8):: Tw, cw !Wall temperature and specific heat
!Neutral Properties
@ -817,16 +811,16 @@ MODULE moduleInput
REAL(8), DIMENSION(:), ALLOCATABLE:: v0
REAL(8):: effTime
REAL(8):: eThreshold !Energy threshold
INTEGER:: speciesID
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
INTEGER:: speciesID, electronSecondaryID
CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, electronSecondary
LOGICAL:: found
INTEGER:: nTypes
CALL config%info('boundary', found, n_children = nBoundary)
ALLOCATE(boundary(1:nBoundary))
DO i = 1, nBoundary
WRITE(istring, '(i2)') i
object = 'boundary(' // TRIM(istring) // ')'
WRITE(iString, '(i2)') i
object = 'boundary(' // TRIM(iString) // ')'
boundary(i)%n = i
CALL config%get(object // '.name', boundary(i)%name, found)
@ -873,8 +867,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)
@ -924,7 +927,6 @@ MODULE moduleInput
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
REAL(8):: volume
CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT
object = 'geometry'
@ -1102,13 +1104,13 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(2):: istring
CHARACTER(2):: iString
INTEGER:: i
CHARACTER(:), ALLOCATABLE:: speciesName
REAL(8), ALLOCATABLE, DIMENSION(:):: r
REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
INTEGER, ALLOCATABLE, DIMENSION(:):: points
REAL(8):: timeStep
REAL(8):: everyTimeStep
CALL config%info('output.probes', found, n_children = nProbes)
@ -1116,7 +1118,7 @@ MODULE moduleInput
DO i = 1, nProbes
WRITE(iString, '(I2)') i
object = 'output.probes(' // trim(istring) // ')'
object = 'output.probes(' // trim(iString) // ')'
CALL config%get(object // '.species', speciesName, found)
CALL config%get(object // '.position', r, found)
@ -1124,16 +1126,14 @@ MODULE moduleInput
CALL config%get(object // '.velocity_2', v2, found)
CALL config%get(object // '.velocity_3', v3, found)
CALL config%get(object // '.points', points, found)
CALL config%get(object // '.timeStep', timeStep, found)
CALL config%get(object // '.timeStep', everyTimeStep, found)
IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)
END DO
CALL resetProbes(tInitial)
END SUBROUTINE readProbes
SUBROUTINE readEMBoundary(config)
@ -1141,7 +1141,6 @@ MODULE moduleInput
USE moduleOutput
USE moduleErrors
USE moduleEM
USE moduleRefParam
USE moduleSpecies
USE json_module
IMPLICIT NONE
@ -1149,34 +1148,72 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(2):: istring
INTEGER:: i, e, s
CHARACTER(:), ALLOCATABLE:: typeEM
REAL(8):: potential
INTEGER:: physicalSurface
CHARACTER(:), ALLOCATABLE:: temporalProfile, temporalProfilePath
INTEGER:: b, s, n, ni
CHARACTER(2):: bString
INTEGER:: info
EXTERNAL:: dgetrf
CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
IF (found) THEN
ALLOCATE(boundaryEM(1:nBoundaryEM))
END IF
DO i = 1, nBoundaryEM
WRITE(istring, '(I2)') i
object = 'boundaryEM(' // trim(istring) // ')'
DO b = 1, nBoundaryEM
WRITE(bString, '(I2)') b
object = 'boundaryEM(' // TRIM(bString) // ')'
CALL config%get(object // '.type', boundEM(i)%typeEM, found)
CALL config%get(object // '.type', typeEM, found)
SELECT CASE(boundEM(i)%typeEM)
SELECT CASE(typeEM)
CASE ("dirichlet")
CALL config%get(object // '.potential', boundEM(i)%potential, found)
IF (.NOT. found) &
CALL config%get(object // '.potential', potential, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
boundEM(i)%potential = boundEM(i)%potential/Volt_ref
CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
IF (.NOT. found) &
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
END IF
CALL config%get(object // '.physicalSurface', physicalSurface, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', &
'readEMBoundary')
END IF
CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential)
CASE ("dirichletTime")
CALL config%get(object // '.potential', potential, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', &
'readEMBoundary')
END IF
CALL config%get(object // '.temporalProfile', temporalProfile, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', &
'readEMBoundary')
END IF
temporalProfilePath = path // temporalProfile
CALL config%get(object // '.physicalSurface', physicalSurface, found)
IF (.NOT. found) THEN
CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', &
'readEMBoundary')
END IF
CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath)
CASE DEFAULT
CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')
END SELECT
@ -1195,18 +1232,28 @@ MODULE moduleInput
END DO
IF (ALLOCATED(boundEM)) THEN
DO e = 1, mesh%numEdges
IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
DO i = 1, nBoundaryEM
IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
CALL boundEM(i)%apply(mesh%edges(e)%obj)
! Modify K matrix due to boundary conditions
DO b = 1, nBoundaryEM
SELECT TYPE(boundary => boundaryEM(b)%obj)
TYPE IS(boundaryEMDirichlet)
DO n = 1, boundary%nNodes
ni = boundary%nodes(n)%obj%n
mesh%K(ni, :) = 0.D0
mesh%K(ni, ni) = 1.D0
END IF
END DO
END IF
END DO
END IF
END DO
TYPE IS(boundaryEMDirichletTime)
DO n = 1, boundary%nNodes
ni = boundary%nodes(n)%obj%n
mesh%K(ni, :) = 0.D0
mesh%K(ni, ni) = 1.D0
END DO
END SELECT
END DO
!Compute the PLU factorization of K once boundary conditions have been read
CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
@ -1227,24 +1274,25 @@ MODULE moduleInput
TYPE(json_file), INTENT(inout):: config
INTEGER:: i
CHARACTER(2):: istring
CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: speciesName
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v
REAL(8), ALLOCATABLE:: T(:), normal(:)
REAL(8), ALLOCATABLE:: temperature(:), normal(:)
REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units
INTEGER:: physicalSurface
INTEGER:: particlesPerEdge
INTEGER:: sp
CALL config%info('inject', found, n_children = nInject)
ALLOCATE(inject(1:nInject))
nPartInj = 0
DO i = 1, nInject
WRITE(istring, '(i2)') i
object = 'inject(' // trim(istring) // ')'
WRITE(iString, '(i2)') i
object = 'inject(' // trim(iString) // ')'
!Find species
CALL config%get(object // '.species', speciesName, found)
@ -1252,7 +1300,7 @@ MODULE moduleInput
CALL config%get(object // '.name', name, found)
CALL config%get(object // '.v', v, found)
CALL config%get(object // '.T', T, found)
CALL config%get(object // '.T', temperature, found)
CALL config%get(object // '.n', normal, found)
IF (.NOT. found) THEN
ALLOCATE(normal(1:3))
@ -1261,8 +1309,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, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
CALL readVelDistr(config, inject(i), object)
@ -1322,28 +1372,28 @@ MODULE moduleInput
TYPE(injectGeneric), INTENT(inout):: inj
CHARACTER(:), ALLOCATABLE, INTENT(in):: object
INTEGER:: i
CHARACTER(2):: istring
CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: fvType
LOGICAL:: found
REAL(8):: v, T, m
REAL(8):: v, temperature, m
!Reads species mass
m = inj%species%m
!Reads distribution functions for velocity
DO i = 1, 3
WRITE(istring, '(i2)') i
CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
WRITE(iString, '(i2)') i
CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found)
IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // &
" found for " // object, 'readVelDistr')
SELECT CASE(fvType)
CASE ("Maxwellian")
T = inj%T(i)
CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
temperature = inj%temperature(i)
CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)
CASE ("Half-Maxwellian")
T = inj%T(i)
CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
temperature = inj%temperature(i)
CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)
CASE ("Delta")
v = inj%vMod*inj%n(i)
@ -1384,5 +1434,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

3
src/modules/mesh/1DCart/moduleMesh1DCart.f90
View file

@ -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
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
!Boundary index

3
src/modules/mesh/1DRad/moduleMesh1DRad.f90
View file

@ -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
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
!Boundary index

16
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

@ -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

53
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -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

3
src/modules/mesh/3DCart/moduleMesh3DCart.f90
View file

@ -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))

21
src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
View file

@ -1,22 +1,22 @@
MODULE moduleMeshOutput0D
CONTAINS
SUBROUTINE printOutput0D(self, t)
SUBROUTINE printOutput0D(self)
USE moduleMesh
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: i
TYPE(outputFormat):: output
CHARACTER(:), ALLOCATABLE:: fileName
DO i = 1, nSpecies
fileName='OUTPUT_' // species(i)%obj%name // '.dat'
IF (t == 0) THEN
IF (timeStep == 0) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", &
"pressure (Pa)", "temperature (K)"
@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj)
WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature
WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density, &
output%velocity, &
output%pressure, &
output%temperature
CLOSE(20)
END DO
END SUBROUTINE printOutput0D
SUBROUTINE printColl0D(self, t)
SUBROUTINE printColl0D(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
CHARACTER(:), ALLOCATABLE:: fileName
INTEGER:: k
fileName='OUTPUT_Collisions.dat'
IF (t == tInitial) THEN
IF (timeStep == tInitial) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions"
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D
END IF
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
CLOSE(20)
END SUBROUTINE printColl0D
SUBROUTINE printEM0D(self, t)
SUBROUTINE printEM0D(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM0D

5
src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
View file

@ -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

38
src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
View file

@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE writeGmsh2FooterElementData
!Prints the scattered properties of particles into the nodes
SUBROUTINE printOutputGmsh2(self, t)
SUBROUTINE printOutputGmsh2(self)
USE moduleMesh
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
time = DBLE(t)*tauMin*ti_ref
time = DBLE(timeStep)*tauMin*ti_ref
DO i = 1, nSpecies
fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t)
fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
CALL writeGmsh2HeaderMesh(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
END DO
CALL writeGmsh2FooterNodeData(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes)
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
END DO
CALL writeGmsh2FooterNodeData(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
END DO
CALL writeGmsh2FooterNodeData(60)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes)
CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO
@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE printOutputGmsh2
!Prints the number of collisions into the volumes
SUBROUTINE printCollGmsh2(self, t)
SUBROUTINE printCollGmsh2(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: numEdges
INTEGER:: k, c
INTEGER:: n
@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2
END SELECT
IF (collOutput) THEN
time = DBLE(t)*tauMin*ti_ref
time = DBLE(timeStep)*tauMin*ti_ref
fileName = formatFileName(prefix, 'Collisions', 'msh', t)
fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2
DO c = 1, interactionMatrix(k)%amount
WRITE(cString, "(I2)") c
title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells)
CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells)
DO n=1, self%numCells
WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c)
END DO
@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2
END SUBROUTINE printCollGmsh2
!Prints the electrostatic EM properties into the nodes and volumes
SUBROUTINE printEMGmsh2(self, t)
SUBROUTINE printEMGmsh2(self)
USE moduleMesh
USE moduleRefParam
USE moduleCaseParam
@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, e
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2
Xi = (/ 0.D0, 0.D0, 0.D0 /)
IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref
time = DBLE(timeStep)*tauMin*ti_ref
fileName = formatFileName(prefix, 'EMField', 'msh', t)
fileName = formatFileName(prefix, 'EMField', 'msh', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (20, file = path // folder // '/' // fileName)
CALL writeGmsh2HeaderMesh(20)
CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes)
CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes)
DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
END DO
CALL writeGmsh2FooterNodeData(20)
CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells)
CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells)
DO e=1, self%numCells
WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref
END DO
CALL writeGmsh2FooterElementData(20)
CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes)
CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes)
DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
END DO

8
src/modules/mesh/inout/moduleMeshInoutCommon.f90
View file

@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon
CHARACTER(LEN=4):: prefix = 'Step'
CONTAINS
PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName)
USE moduleOutput
IMPLICIT NONE
CHARACTER(*), INTENT(in):: prefix, suffix, extension
INTEGER, INTENT(in), OPTIONAL:: t
INTEGER, INTENT(in), OPTIONAL:: timeStep
CHARACTER (LEN=iterationDigits):: tString
CHARACTER(:), ALLOCATABLE:: fileName
IF (PRESENT(t)) THEN
WRITE(tString, iterationFormat) t
IF (PRESENT(timeStep)) THEN
WRITE(tString, iterationFormat) timeStep
fileName = prefix // '_' // tString // '_' // suffix // '.' // extension
ELSE

5
src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
View file

@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU
CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r(1:3) !3 generic coordinates
INTEGER:: fileID, error, found
INTEGER:: fileID
CHARACTER(LEN=256):: line
INTEGER:: numNodes, numElements, numEdges
INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
@ -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

46
src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
View file

@ -209,23 +209,22 @@ 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>'
END SUBROUTINE writeEM
SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection)
USE moduleCaseParam
USE moduleOutput
USE moduleRefParam
IMPLICIT NONE
INTEGER:: fileID
INTEGER, INTENT(in):: t
CHARACTER(*):: fileNameStep, fileNameCollection
IF (t == tInitial) THEN
IF (timeStep == tInitial) THEN
!Create collection file
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection)
@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU
!Write iteration file in collection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND')
WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") &
'<DataSet timestep="', DBLE(timeStep)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
!Close collection file
IF (t == tFinal) THEN
IF (timeStep == tFinal) THEN
WRITE (fileID + 1, "(2X, A)") '</Collection>'
WRITE (fileID + 1, "(A)") '</VTKFile>'
@ -307,22 +307,21 @@ MODULE moduleMeshOutputVTU
END SUBROUTINE writeAverage
SUBROUTINE printOutputVTU(self,t)
SUBROUTINE printOutputVTU(self)
USE moduleMesh
USE moduleSpecies
USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
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
DO i = 1, nSpecies
fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
@ -338,29 +337,27 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
CALL writeCollection(fileID, fileName, filenameCollection)
END DO
END SUBROUTINE printOutputVTU
SUBROUTINE printCollVTU(self,t)
SUBROUTINE printCollVTU(self)
USE moduleMesh
USE moduleOutput
USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
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
IF (collOutput) THEN
fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
@ -376,26 +373,26 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
CALL writeCollection(fileID, fileName, filenameCollection)
END IF
END SUBROUTINE printCollVTU
SUBROUTINE printEMVTU(self, t)
SUBROUTINE printEMVTU(self)
USE moduleMesh
USE moduleMeshInoutCommon
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 64
IF (emOutput) THEN
fileName = formatFileName(prefix, 'EMField', 'vtu', t)
fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
@ -411,7 +408,7 @@ MODULE moduleMeshOutputVTU
!Write collection file for time plotting
fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
CALL writeCollection(fileID, fileName, filenameCollection)
END IF
@ -424,9 +421,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

40
src/modules/mesh/moduleMesh.f90
View file

@ -59,6 +59,13 @@ MODULE moduleMesh
END TYPE meshNodeCont
! Array of pointers to nodes.
TYPE:: meshNodePointer
CLASS(meshNode), POINTER:: obj
CONTAINS
END TYPE meshNodePointer
!Type for array of boundary functions (one per species)
TYPE, PUBLIC:: fBoundaryGeneric
PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
@ -76,8 +83,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
@ -372,10 +379,9 @@ MODULE moduleMesh
END SUBROUTINE connectMesh_interface
!Prints number of collisions in each cell
SUBROUTINE printColl_interface(self, t)
SUBROUTINE printColl_interface(self)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface
@ -403,18 +409,16 @@ MODULE moduleMesh
ABSTRACT INTERFACE
!Prints Species data
SUBROUTINE printOutput_interface(self, t)
SUBROUTINE printOutput_interface(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface
!Prints EM info
SUBROUTINE printEM_interface(self, t)
SUBROUTINE printEM_interface(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface
@ -613,6 +617,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 +628,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
@ -787,7 +793,7 @@ MODULE moduleMesh
END FUNCTION findCellBrute
!Computes collisions in element
SUBROUTINE doCollisions(self, t)
SUBROUTINE doCollisions(self)
USE moduleCollisions
USE moduleSpecies
USE moduleList
@ -795,10 +801,10 @@ MODULE moduleMesh
USE moduleRandom
USE moduleOutput
USE moduleMath
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self
INTEGER, INTENT(in):: t
INTEGER:: e
CLASS(meshCell), POINTER:: cell
INTEGER:: k, i, j
@ -814,7 +820,7 @@ MODULE moduleMesh
REAL(8):: rnd_real !Random number for collision
INTEGER:: rnd_int !Random number for collision
IF (MOD(t, everyColl) == 0) THEN
IF (MOD(timeStep, everyColl) == 0) THEN
!Collisions need to be performed in this iteration
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
DO e=1, self%numCells
@ -1023,6 +1029,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
@ -1107,6 +1116,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

8
src/modules/mesh/moduleMeshBoundary.f90
View file

@ -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))

24
src/modules/moduleBoundary.f90
View file

@ -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
@ -103,17 +104,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
@ -126,7 +129,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)

223
src/modules/moduleInject.f90
View file

@ -54,15 +54,16 @@ MODULE moduleInject
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: vMod !Velocity (module)
REAL(8):: T(1:3) !Temperature
REAL(8):: temperature(1:3) !Temperature
REAL(8):: n(1:3) !Direction of injection
LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
INTEGER:: nParticles !Number of particles to introduce each time step
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, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
USE moduleMesh
USE moduleRefParam
USE moduleConstParam
@ -86,50 +87,29 @@ 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
REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
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
! TODO: Make this only available for charge species
self%nParticles = INT(flow*tauInject*ti_ref/(qe*abs(species(sp)%obj%qm*species(sp)%obj%m)*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%temperature = temperature / 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
@ -161,15 +141,78 @@ 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
SELECT TYPE(sp => self%species)
CLASS IS(speciesCharged)
fluxPerStep = flow/(qe*abs(sp%q))
CLASS DEFAULT
call criticalError('Attempted to assign a flux in "A" to a species without charge.', 'initInject')
END SELECT
CASE ("Am2")
!Input current in Ampers per square meter
SELECT TYPE(sp => self%species)
CLASS IS(speciesCharged)
fluxPerStep = flow*self%surface*L_ref**2/(qe*abs(sp%q))
CLASS DEFAULT
call criticalError('Attempted to assign a flux in "Am2" to a species without charge.', 'initInject')
END SELECT
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
@ -204,23 +247,23 @@ MODULE moduleInject
END SUBROUTINE doInjects
SUBROUTINE initVelDistMaxwellian(velDist, T, m)
SUBROUTINE initVelDistMaxwellian(velDist, temperature, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: T, m
REAL(8), INTENT(in):: temperature, m
velDist = velDistMaxwellian(vTh = DSQRT(T/m))
velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: T, m
REAL(8), INTENT(in):: temperature, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m))
velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian
@ -283,9 +326,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, e
INTEGER:: n, sp
CLASS(meshEdge), POINTER:: randomEdge
REAL(8):: direction(1:3)
@ -300,59 +342,62 @@ 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
partInj(n)%v = 0.D0
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /)
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /)
!If velocity is not in the right direction, invert it
IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
partInj(n)%v = - partInj(n)%v
!If velocity is not in the right direction, invert it
IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
partInj(n)%v = - partInj(n)%v
END IF
END IF
!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

54
src/modules/moduleProbe.f90
View file

@ -27,7 +27,7 @@ MODULE moduleProbe
CONTAINS
!Functions for probeDistFunc type
SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep)
SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep)
USE moduleCaseParam
USE moduleRefParam
USE moduleSpecies
@ -41,7 +41,7 @@ MODULE moduleProbe
REAL(8), INTENT(in):: r(1:3)
REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
INTEGER, INTENT(in):: points(1:3)
REAL(8), INTENT(in):: timeStep
REAL(8), INTENT(in):: everyTimeStep
INTEGER:: sp, i
REAL(8):: dv(1:3)
@ -91,17 +91,17 @@ MODULE moduleProbe
1:self%nv(3)))
!Number of iterations between output
IF (timeStep == 0.D0) THEN
IF (everyTimeStep == 0.D0) THEN
self%every = 1
ELSE
self%every = NINT(timeStep/ tauMin / ti_ref)
self%every = NINT(everyTimeStep/ tauMin / ti_ref)
END IF
!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,40 +162,40 @@ 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
END SUBROUTINE calculate
SUBROUTINE output(self, t)
SUBROUTINE output(self)
USE moduleOutput
USE moduleRefParam
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
CLASS(probeDistFunc), INTENT(inout):: self
INTEGER, INTENT(in):: t
CHARACTER (LEN=iterationDigits):: tstring
CHARACTER (LEN=3):: pstring
CHARACTER(:), ALLOCATABLE:: filename
@ -204,14 +204,14 @@ MODULE moduleProbe
!Divide by the velocity cube volume
self%f = self%f * self%dvInv
WRITE(tstring, iterationFormat) t
WRITE(tstring, iterationFormat) timeStep
WRITE(pstring, "(I3.3)") self%id
fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (10, file = path // folder // '/' // fileName)
WRITE(10, "(A1, 1X, A)") "# ", self%species%name
WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(t)*tauMin*ti_ref, " s"
WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(timeStep)*tauMin*ti_ref, " s"
WRITE(10, "(A1, A19, 3(A20))") "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
DO i = 1, self%nv(1)
DO j = 1, self%nv(2)
@ -252,15 +252,15 @@ MODULE moduleProbe
END SUBROUTINE doProbes
SUBROUTINE outputProbes(t)
SUBROUTINE outputProbes()
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: i
DO i = 1, nProbes
IF (probe(i)%update) THEN
CALL probe(i)%output(t)
CALL probe(i)%output()
END IF
@ -268,15 +268,15 @@ MODULE moduleProbe
END SUBROUTINE outputProbes
SUBROUTINE resetProbes(t)
SUBROUTINE resetProbes()
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: i
DO i = 1, nProbes
probe(i)%f = 0.D0
probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0
probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0
END DO

6
src/modules/output/moduleOutput.f90
View file

@ -160,12 +160,12 @@ MODULE moduleOutput
END SUBROUTINE calculateOutput
SUBROUTINE printTime(t, first)
SUBROUTINE printTime(first)
USE moduleSpecies
USE moduleCompTime
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
INTEGER, INTENT(in):: t
LOGICAL, INTENT(in), OPTIONAL:: first
CHARACTER(:), ALLOCATABLE:: fileName
@ -187,7 +187,7 @@ MODULE moduleOutput
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
CLOSE(20)

200
src/modules/solver/electromagnetic/moduleEM.f90
View file

@ -1,52 +1,200 @@
!Module to solve the electromagnetic field
MODULE moduleEM
USE moduleMesh
USE moduleTable
IMPLICIT NONE
TYPE:: boundaryEM
CHARACTER(:), ALLOCATABLE:: typeEM
INTEGER:: physicalSurface
! Generic type for electromagnetic boundary conditions
TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric
INTEGER:: nNodes
TYPE(meshNodePointer), ALLOCATABLE:: nodes(:)
CONTAINS
PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply
!PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better.
END TYPE boundaryEMGeneric
ABSTRACT INTERFACE
! Apply boundary condition to the load vector for the Poission equation
SUBROUTINE applyEM_interface(self, vectorF)
IMPORT boundaryEMGeneric
CLASS(boundaryEMGeneric), INTENT(in):: self
REAL(8), INTENT(inout):: vectorF(:)
END SUBROUTINE applyEM_interface
END INTERFACE
TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet
REAL(8):: potential
CONTAINS
PROCEDURE, PASS:: apply
! boundaryEMGeneric DEFERRED PROCEDURES
PROCEDURE, PASS:: apply => applyDirichlet
END TYPE boundaryEM
END TYPE boundaryEMDirichlet
TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
REAL(8):: potential
TYPE(table1D):: temporalProfile
CONTAINS
! boundaryEMGeneric DEFERRED PROCEDURES
PROCEDURE, PASS:: apply => applyDirichletTime
END TYPE boundaryEMDirichletTime
! Container for boundary conditions
TYPE:: boundaryEMCont
CLASS(boundaryEMGeneric), ALLOCATABLE:: obj
END TYPE boundaryEMCont
INTEGER:: nBoundaryEM
TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:)
!Information of charge and reference parameters for rho vector
REAL(8), ALLOCATABLE:: qSpecies(:)
CONTAINS
!Apply boundary conditions to the K matrix for Poisson's equation
SUBROUTINE apply(self, edge)
SUBROUTINE findNodes(self, physicalSurface)
USE moduleMesh
IMPLICIT NONE
CLASS(boundaryEM), INTENT(in):: self
CLASS(meshEdge):: edge
INTEGER:: nNodes
INTEGER, ALLOCATABLE:: nodes(:)
INTEGER:: n
CLASS(boundaryEMGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: physicalSurface
CLASS(meshEdge), POINTER:: edge
INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
INTEGER:: nNodes, nodesNew
INTEGER:: e, n
nNodes = edge%nNodes
nodes = edge%getNodes(nNodes)
!Temporal array to hold nodes
ALLOCATE(nodes(0))
DO n = 1, nNodes
SELECT CASE(self%typeEM)
CASE ("dirichlet")
mesh%K(nodes(n), :) = 0.D0
mesh%K(nodes(n), nodes(n)) = 1.D0
mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
mesh%nodes(nodes(n))%obj%emData%phi = self%potential
! Loop thorugh the edges and identify those that are part of the boundary
DO e = 1, mesh%numEdges
edge => mesh%edges(e)%obj
IF (edge%physicalSurface == physicalSurface) THEN
! Edge is of the right boundary index
! Get nodes in the edge
nNodes = edge%nNodes
nodesEdge = edge%getNodes(nNodes)
! Collect all nodes that are not already in the temporal array
DO n = 1, nNodes
IF (ANY(nodes == nodesEdge(n))) THEN
! Node already in array, skip
CYCLE
END SELECT
ELSE
! If not, add element to array of nodes
nodes = [nodes, nodesEdge(n)]
END IF
END DO
END IF
END DO
END SUBROUTINE apply
! Point boundary to nodes
nNodes = SIZE(nodes)
ALLOCATE(self%nodes(nNodes))
self%nNodes = nNodes
DO n = 1, nNodes
self%nodes(n)%obj => mesh%nodes(nodes(n))%obj
END DO
END SUBROUTINE findNodes
! Initialize Dirichlet boundary condition
SUBROUTINE initDirichlet(self, physicalSurface, potential)
USE moduleRefParam, ONLY: Volt_ref
IMPLICIT NONE
CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
INTEGER, INTENT(in):: physicalSurface
REAL(8), INTENT(in):: potential
! Allocate boundary edge
ALLOCATE(boundaryEMDirichlet:: self)
SELECT TYPE(self)
TYPE IS(boundaryEMDirichlet)
self%potential = potential / Volt_ref
CALL findNodes(self, physicalSurface)
END SELECT
END SUBROUTINE initDirichlet
! Initialize Dirichlet boundary condition
SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile)
USE moduleRefParam, ONLY: Volt_ref, ti_ref
IMPLICIT NONE
CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self
INTEGER, INTENT(in):: physicalSurface
REAL(8), INTENT(in):: potential
CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile
! Allocate boundary edge
ALLOCATE(boundaryEMDirichletTime:: self)
SELECT TYPE(self)
TYPE IS(boundaryEMDirichletTime)
self%potential = potential / Volt_ref
CALL findNodes(self, physicalSurface)
CALL self%temporalProfile%init(temporalProfile)
CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
END SELECT
END SUBROUTINE initDirichletTime
!Apply Dirichlet boundary condition to the poisson equation
SUBROUTINE applyDirichlet(self, vectorF)
USE moduleMesh
IMPLICIT NONE
CLASS(boundaryEMDirichlet), INTENT(in):: self
REAL(8), INTENT(inout):: vectorF(:)
INTEGER:: n, ni
DO n = 1, self%nNodes
self%nodes(n)%obj%emData%phi = self%potential
vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
END DO
END SUBROUTINE applyDirichlet
!Apply Dirichlet boundary condition with time temporal profile
SUBROUTINE applyDirichletTime(self, vectorF)
USE moduleMesh
USE moduleCaseParam, ONLY: timeStep, tauMin
IMPLICIT NONE
CLASS(boundaryEMDirichletTime), INTENT(in):: self
REAL(8), INTENT(inout):: vectorF(:)
REAL(8):: timeFactor
INTEGER:: n, ni
timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
DO n = 1, self%nNodes
self%nodes(n)%obj%emData%phi = self%potential * timeFactor
vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
END DO
END SUBROUTINE applyDirichletTime
!Assemble the source vector based on the charge density to solve Poisson's equation
SUBROUTINE assembleSourceVector(vectorF, n_e)
@ -60,7 +208,7 @@ MODULE moduleEM
REAL(8), ALLOCATABLE:: rho(:)
REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes)
INTEGER:: nNodes
INTEGER:: e, i, ni
INTEGER:: e, i, ni, b
CLASS(meshNode), POINTER:: node
! !$OMP SINGLE

42
src/modules/solver/moduleSolver.f90
View file

@ -493,52 +493,46 @@ MODULE moduleSolver
END SUBROUTINE updateParticleCell
!Update the information about if a species needs to be moved this iteration
SUBROUTINE updatePushSpecies(self, t)
SUBROUTINE updatePushSpecies(self)
USE moduleSpecies
USE moduleCaseparam, ONLY: timeStep
IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t
INTEGER:: s
DO s=1, nSpecies
self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0
self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0
END DO
END SUBROUTINE updatePushSpecies
!Output the different data and information
SUBROUTINE doOutput(t)
SUBROUTINE doOutput()
USE moduleMesh
USE moduleOutput
USE moduleSpecies
USE moduleCompTime
USE moduleProbe
USE moduleCaseParam, ONLY: timeStep
IMPLICIT NONE
INTEGER, INTENT(in):: t
IF (t == tInitial) THEN
CALL SYSTEM('git rev-parse HEAD > ' // path // folder // '/' // 'fpack_commit.txt')
END IF
CALL outputProbes(t)
CALL outputProbes()
counterOutput = counterOutput + 1
IF (counterOutput >= triggerOutput .OR. &
t == tFinal .OR. t == tInitial) THEN
timeStep == tFinal .OR. timeStep == tInitial) THEN
!Resets output counter
counterOutput=0
CALL mesh%printOutput(t)
IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
CALL mesh%printEM(t)
WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
CALL mesh%printOutput()
IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl()
CALL mesh%printEM()
WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal
WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld
IF (t == 0) THEN
IF (timeStep == 0) THEN
WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms"
ELSE
@ -556,34 +550,32 @@ MODULE moduleSolver
counterCPUTime = counterCPUTime + 1
IF (counterCPUTime >= triggerCPUTime .OR. &
t == tFinal .OR. t == tInitial) THEN
timeStep == tFinal .OR. timeStep == tInitial) THEN
!Reset CPU Time counter
counterCPUTime = 0
CALL printTime(t, t == 0)
CALL printTime(timeStep == 0)
END IF
!Output average values
IF (useAverage .AND. t == tFinal) THEN
IF (useAverage .AND. timeStep == tFinal) THEN
CALL mesh%printAverage()
END IF
END SUBROUTINE doOutput
SUBROUTINE doAverage(t)
SUBROUTINE doAverage()
USE moduleAverage
USE moduleMesh
IMPLICIT NONE
INTEGER, INTENT(in):: t
INTEGER:: tAverage, n
IF (useAverage) THEN
tAverage = t - tAverageStart
tAverage = timeStep - tAverageStart
IF (tAverage == 1) THEN
!First iteration in which average scheme is used