diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index c3b66ab..f8ab920 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index 1b4e258..75f0106 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -401,6 +401,9 @@ make \item \textbf{bounraiesParticle}: Logical. Determines if the properties of the boundaries are written. \textit{Warning}: Not all particle models produce this output, only those whose properties change during the simulation. + \item \textbf{injects}: Logical. + Determines if the properties of the injection of particles are written. + \textit{Warning}: Not all injection types produce this output, only those whose properties change during the simulation. \item \textbf{probes}: Array of objects. Defines the probes employed for obtaining the distribution function at specific positions. See Sec.~\ref{sec:probing} for more information. @@ -688,6 +691,14 @@ make Units of $\unit{K}$ Array dimension $3$. Temperature in each direction. + \item \textbf{type}: Character. + Type of injection. + Available values are: + \begin{itemize} + \item \textbf{constant}: The flow does not change. + This is the default value if none is provided. + \item \textbf{quasiNeutral}: The weight of the injected species' particles is changed over time to keep a zero electric field at the boundary. + \end{itemize} \item \textbf{physicalSurface}: Integer. Identification of the edge in the mesh file. \item \textbf{particlesPerEdge}: Integer. diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 921d4ac..cb68754 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -1,15 +1,39 @@ ! FPAKC main program PROGRAM fpakc - USE moduleCompTime - USE moduleCaseParam - USE moduleInput - USE moduleInject - USE moduleSolver - USE moduleMesh - USE moduleProbe - USE moduleErrors - USE OMP_LIB - IMPLICIT NONE + use moduleCompTime, only: tStep, & + tEMField, & + tCoul, & + tColl, & + tPush, & + tReset, & + tWeight + use omp_lib, only: omp_get_wtime + use moduleErrors, only: criticalError, & + verboseError + use moduleInput, only: readCOnfig, & + initOutput + use moduleCaseParam, only: timeStep, & + tInitial, & + tFinal + use moduleProbe, only: resetProbes + use moduleSolver, only: doScatter, & + doEMField, & + doOutput, & + solver, & + doPushes, & + doReset, & + doAverage, & + doInjects + use moduleMesh, only: boundariesEM_update, & + boundariesEM_update, & + boundariesParticle_update, & + mesh, & + meshForMCC, & + doMCCollisions, & + doCollisions, & + doCoulombScattering + use moduleInject, only: updateInjects + implicit none ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 @@ -71,6 +95,9 @@ PROGRAM fpakc ! Update EM boundary models call boundariesEM_update() + ! Update injects + call updateInjects() + !Checks if a species needs to me moved in this iteration CALL solver%updatePushSpecies() diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 1b8ec67..a1bc3fc 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -53,7 +53,7 @@ MODULE moduleRandom implicit none real(8):: rnd - real(8):: v1, v2, Rsquare + real(8):: v1, v2 v1 = 0.d0 do while (v1 <= 0.d0) @@ -67,19 +67,20 @@ MODULE moduleRandom end function randomMaxwellian !Returns a random number in a Maxwellian distribution of mean 0 and width 1 - FUNCTION randomHalfMaxwellian() RESULT(rnd) + function randomHalfMaxwellian() result(rnd) IMPLICIT NONE - REAL(8):: rnd - REAL(8):: x + real(8):: rnd + real(8):: v1 rnd = 0.D0 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - END DO + v1 = 0.D0 + do while (v1 == 0.D0) + v1 = random() + + end do - rnd = DSQRT(-DLOG(x)) + rnd = sqrt(-2.d0*log(v1)) END FUNCTION randomHalfMaxwellian diff --git a/src/modules/common/velocityDistribution.f90 b/src/modules/common/velocityDistribution.f90 index 5b6c342..359eebd 100644 --- a/src/modules/common/velocityDistribution.f90 +++ b/src/modules/common/velocityDistribution.f90 @@ -53,12 +53,11 @@ module velocityDistribution function randomVelMaxwellian(self) result (v) use moduleRandom implicit none - class(velDistMaxwellian), intent(in):: self real(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian()/sqrt(2.d0) + v = self%vTh*randomMaxwellian() end function randomVelMaxwellian diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 8b7aaff..eb8bde5 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -493,6 +493,7 @@ MODULE moduleInput CALL config%get(object // '.numColl', collOutput, found) CALL config%get(object // '.EMField', emOutput, found) call config%get(object // '.boundariesParticle', boundaryParticleOutput, found) + call config%get(object // '.injects', injectOutput, found) CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found) IF (.NOT. found) THEN @@ -628,6 +629,7 @@ MODULE moduleInput USE moduleMesh USE moduleCaseParam USE moduleRefParam + use moduleOutput, only: collOutput USE OMP_LIB USE json_module IMPLICIT NONE @@ -1426,53 +1428,19 @@ MODULE moduleInput CHARACTER(2):: iString CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: speciesName - CHARACTER(:), ALLOCATABLE:: name - REAL(8):: v - 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 + ALLOCATE(injects(1:nInject)) DO i = 1, nInject WRITE(iString, '(i2)') i object = 'inject(' // trim(iString) // ')' - !Find species - CALL config%get(object // '.species', speciesName, found) - sp = speciesName2Index(speciesName) + CALL initInject(injects(i)%obj, i, object, config) - CALL config%get(object // '.name', name, found) - CALL config%get(object // '.v', v, found) - CALL config%get(object // '.T', temperature, found) - CALL config%get(object // '.n', normal, found) - IF (.NOT. found) THEN - ALLOCATE(normal(1:3)) - normal = 0.D0 - END IF - 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, temperature, flow, units, sp, physicalSurface, particlesPerEdge) - - CALL readVelDistr(config, inject(i), object) + CALL readVelDistr(config, injects(i)%obj, object) END DO - !Allocate array for injected particles - IF (nPartInj > 0) THEN - ALLOCATE(partInj(1:nPartInj)) - - END IF - END SUBROUTINE readInject SUBROUTINE readAverage(config) @@ -1522,7 +1490,7 @@ MODULE moduleInput IMPLICIT NONE TYPE(json_file), INTENT(inout):: config - TYPE(injectGeneric), INTENT(inout):: inj + class(injectGeneric), INTENT(inout):: inj CHARACTER(:), ALLOCATABLE, INTENT(in):: object INTEGER:: i CHARACTER(2):: iString diff --git a/src/modules/makefile b/src/modules/makefile index 9f62c85..5789276 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -14,7 +14,7 @@ output.o: moduleSpecies.o common.o mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o $(MAKE) -C mesh all -solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o +solver.o: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o $(MAKE) -C solver all init.o: common.o solver.o moduleInject.o diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 0ca61ad..acd4189 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -1,4 +1,5 @@ MODULE moduleMeshOutputGmsh2 + use moduleOutput, only: prefix, informFileCreation, formatFileName, generateFilePath CONTAINS !Header for mesh format @@ -16,7 +17,7 @@ MODULE moduleMeshOutputGmsh2 !Node data subroutines !Header SUBROUTINE writeGmsh2HeaderNodeData(fileID, title, iteration, time, dimensions, nNodes) - use moduleOutput + use moduleOutput, only: fmtInt, fmtReal IMPLICIT NONE INTEGER, INTENT(in):: fileID @@ -50,7 +51,7 @@ MODULE moduleMeshOutputGmsh2 !Element data subroutines !Header SUBROUTINE writeGmsh2HeaderElementData(fileID, title, iteration, time, dimensions, nVols) - use moduleOutput + use moduleOutput, only: fmtInt, fmtReal IMPLICIT NONE INTEGER, INTENT(in):: fileID @@ -86,7 +87,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleMesh USE moduleRefParam USE moduleSpecies - USE moduleOutput + use moduleOutput, only: outputFormat, calculateOutput, fmtReal USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE @@ -141,7 +142,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleRefParam USE moduleCaseParam USE moduleCollisions - USE moduleOutput + use moduleOutput, only: collOutput, fmtInt IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self @@ -199,7 +200,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleMesh USE moduleRefParam USE moduleCaseParam - USE moduleOutput + use moduleOutput, only: emOutput IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -248,7 +249,6 @@ MODULE moduleMeshOutputGmsh2 USE moduleMesh USE moduleRefParam USE moduleSpecies - USE moduleOutput USE moduleAverage IMPLICIT NONE diff --git a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 index 7f95788..74ac0ea 100644 --- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 @@ -1,9 +1,14 @@ module moduleMeshOutputText + use moduleOutput, only: prefix, & + informFileCreation, & + formatFileName, & + generateFilePath + contains subroutine writeSpeciesOutput(self, fileID, speciesIndex) use moduleMesh - use moduleOutput + use moduleOutput, only: calculateOutput, outputFormat use moduleRefParam, only: L_ref implicit none @@ -88,7 +93,6 @@ module moduleMeshOutputText speciesIndex) use moduleMesh - use moduleOutput use moduleAverage use moduleRefParam, only: L_ref implicit none @@ -151,8 +155,8 @@ module moduleMeshOutputText subroutine printCollText(self) use moduleMesh - use moduleOutput use moduleCaseParam, only: timeStep + use moduleOutput, only: collOutput implicit none class(meshGeneric), intent(in):: self @@ -190,6 +194,7 @@ module moduleMeshOutputText subroutine printEMText(self) use moduleMesh use moduleCaseParam, only: timeStep + use moduleOutput, only: emOutput implicit none class(meshParticles), intent(in):: self @@ -236,7 +241,7 @@ module moduleMeshOutputText write(fileIDMean, '(5(A,","),A)') '"Position (m)"', & '"Density, mean (m^-3)"', & - '"Velocity, mean (m s^-1):0"', '"Velocity (m s^-1):1"', '"Velocity (m s^-1):2"', & + '"Velocity, mean (m s^-1):0"', '"Velocity, mean (m s^-1):1"', '"Velocity, mean (m s^-1):2"', & '"Temperature, mean (K)"' fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv') @@ -245,7 +250,7 @@ module moduleMeshOutputText write(fileIDDeviation, '(5(A,","),A)') '"Position (m)"', & '"Density, deviation (m^-3)"', & - '"Velocity, deviation (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (m s^-1):2"', & + '"Velocity, deviation (m s^-1):0"', 'Velocity, deviation (m s^-1):1"', 'Velocity, deviation (m s^-1):2"', & '"Temperature, deviation (K)"' call writeAverage(self, fileIDMean, fileIDDeviation, s) diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 index 3462940..8058896 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -1,9 +1,11 @@ MODULE moduleMeshOutputVTU + use moduleOutput, only: prefix, informFileCreation, formatFileName, generateFilePath CONTAINS SUBROUTINE writeHeader(nNodes, nCells, fileID) USE moduleMesh + use moduleOutput, only: fmtInt IMPLICIT NONE INTEGER, INTENT(in):: nNodes, nCells @@ -77,6 +79,7 @@ MODULE moduleMeshOutputVTU SUBROUTINE writeMesh(self, fileID) USE moduleMesh USE moduleRefParam + use moduleOutput, only: fmtReal, fmtInt IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self @@ -159,6 +162,7 @@ MODULE moduleMeshOutputVTU SUBROUTINE writeCollOutput(self,fileID) USE moduleMesh USE moduleCollisions + use moduleOutput, only: fmtInt IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self @@ -185,6 +189,7 @@ MODULE moduleMeshOutputVTU SUBROUTINE writeEM(self, fileID) USE moduleMesh USE moduleRefParam + use moduleOutput, only: fmtReal IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -381,6 +386,7 @@ MODULE moduleMeshOutputVTU SUBROUTINE printEMVTU(self) USE moduleMesh USE moduleCaseParam, ONLY: timeStep + use moduleOutput, only: emOutput IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -414,7 +420,6 @@ MODULE moduleMeshOutputVTU SUBROUTINE printAverageVTU(self) USE moduleMesh - use moduleOutput USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 7ce7565..6cc48ee 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -1,7 +1,7 @@ !moduleMesh: General module for Finite Element mesh MODULE moduleMesh USE moduleList - USE moduleOutput + USE moduleOutput, only: outputNode, emNode USE moduleCollisions use moduleSpecies, only: nSpecies @@ -663,50 +663,6 @@ MODULE moduleMesh end type boundaryParticleGeneric - interface - ! Init interfaces - ! wallTemeprature - module subroutine wallTemperature_init(boundary, T, c) - CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary - REAL(8), INTENT(in):: T, c !Wall temperature and specific heat - - end subroutine wallTemperature_init - - module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) - class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - real(8), intent(in):: mImpact - real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties - integer, intent(in):: ion - real(8), intent(in):: effTime - character(:), allocatable, intent(in):: crossSection - real(8), intent(in):: eThreshold - integer, optional, intent(in):: electronSecondary - - end subroutine ionization_init - - ! quasiNeutrality - module subroutine quasiNeutrality_init(boundary, s_incident) - class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - integer, intent(in):: s_incident - - end subroutine quasiNeutrality_init - - ! outflowAdaptive - module subroutine outflowAdaptive_init(boundary, s_incident) - class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - integer, intent(in):: s_incident - - end subroutine outflowAdaptive_init - - ! Get the index of the species from its name - module function boundaryParticleName_to_Index(boundaryName) result(bp) - character(:), allocatable:: boundaryName - integer:: bp - - end function boundaryParticleName_to_Index - - end interface - abstract interface ! Apply the effects of the boundary to the particle subroutine applyParticle_interface(self, edge, part) @@ -745,6 +701,18 @@ MODULE moduleMesh END TYPE boundaryReflection + interface + module subroutine reflection_apply(self, edge, part) + use moduleSpecies + + class(boundaryReflection), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine reflection_apply + + end interface + !Absorption boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAbsorption CONTAINS @@ -752,6 +720,18 @@ MODULE moduleMesh END TYPE boundaryAbsorption + interface + module subroutine absorption_apply(self, edge, part) + use moduleSpecies + + class(boundaryAbsorption), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine absorption_apply + + end interface + !Transparent boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryTransparent CONTAINS @@ -759,6 +739,18 @@ MODULE moduleMesh END TYPE boundaryTransparent + interface + module subroutine transparent_apply(self, edge, part) + use moduleSpecies + + class(boundaryTransparent), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine transparent_apply + + end interface + !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAxis CONTAINS @@ -766,6 +758,18 @@ MODULE moduleMesh END TYPE boundaryAxis + interface + module subroutine symmetryAxis_apply(self, edge, part) + use moduleSpecies + + class(boundaryAxis), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine symmetryAxis_apply + + end interface + !Wall Temperature boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryWallTemperature !Thermal velocity of the wall: square root(Wall temperature X specific heat) @@ -775,6 +779,24 @@ MODULE moduleMesh END TYPE boundaryWallTemperature + interface + module subroutine wallTemperature_init(boundary, T, c) + CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary + REAL(8), INTENT(in):: T, c !Wall temperature and specific heat + + end subroutine wallTemperature_init + + module subroutine wallTemperature_apply(self, edge, part) + use moduleSpecies + + class(boundaryWallTemperature), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine wallTemperature_apply + + end interface + !Ionization boundary TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryIonization REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. @@ -791,6 +813,30 @@ MODULE moduleMesh END TYPE boundaryIonization + interface + module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + real(8), intent(in):: mImpact + real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties + integer, intent(in):: ion + real(8), intent(in):: effTime + character(:), allocatable, intent(in):: crossSection + real(8), intent(in):: eThreshold + integer, optional, intent(in):: electronSecondary + + end subroutine ionization_init + + module subroutine ionization_apply(self, edge, part) + use moduleSpecies + + class(boundaryIonization), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine ionization_apply + + end interface + ! Ensures quasi-neutrality by changing the reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality real(8), allocatable:: alpha(:) ! Reflection parameter @@ -801,6 +847,24 @@ MODULE moduleMesh end type boundaryQuasiNeutrality + interface + module subroutine quasiNeutrality_init(boundary, s_incident) + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident + + end subroutine quasiNeutrality_init + + module subroutine quasiNeutrality_apply(self, edge, part) + use moduleSpecies + + class(boundaryQuasiNeutrality), intent(inout):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine quasiNeutrality_apply + + end interface + !Boundary for quasi-neutral outflow adjusting reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive real(8), allocatable:: velocity_shift(:) @@ -812,68 +876,11 @@ MODULE moduleMesh end type boundaryOutflowAdaptive interface - module subroutine reflection_apply(self, edge, part) - use moduleSpecies + module subroutine outflowAdaptive_init(boundary, s_incident) + class(boundaryParticleGeneric), allocatable, intent(inout):: boundary + integer, intent(in):: s_incident - class(boundaryReflection), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine reflection_apply - - module subroutine absorption_apply(self, edge, part) - use moduleSpecies - - class(boundaryAbsorption), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine absorption_apply - - module subroutine transparent_apply(self, edge, part) - use moduleSpecies - - class(boundaryTransparent), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine transparent_apply - - module subroutine symmetryAxis_apply(self, edge, part) - use moduleSpecies - - class(boundaryAxis), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine symmetryAxis_apply - - module subroutine wallTemperature_apply(self, edge, part) - use moduleSpecies - - class(boundaryWallTemperature), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine wallTemperature_apply - - module subroutine ionization_apply(self, edge, part) - use moduleSpecies - - class(boundaryIonization), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine ionization_apply - - module subroutine quasiNeutrality_apply(self, edge, part) - use moduleSpecies - - class(boundaryQuasiNeutrality), intent(inout):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - - end subroutine quasiNeutrality_apply + end subroutine outflowAdaptive_init module subroutine outflowAdaptive_apply(self, edge, part) use moduleSpecies @@ -884,7 +891,10 @@ MODULE moduleMesh end subroutine outflowAdaptive_apply - ! Generic basic boundary conditions to use internally in the code + end interface + + ! Generic basic boundary conditions to use internally in the code + interface module subroutine genericReflection(edge, part) use moduleSpecies @@ -927,6 +937,13 @@ MODULE moduleMesh type(boundaryParticleCont), allocatable, target:: boundariesParticle(:) interface + ! Get the index of the species from its name + module function boundaryParticleName_to_Index(boundaryName) result(bp) + character(:), allocatable:: boundaryName + integer:: bp + + end function boundaryParticleName_to_Index + ! Update the particle boundary models module subroutine boundariesParticle_update() @@ -952,33 +969,6 @@ MODULE moduleMesh end type boundaryEMGeneric - interface - module subroutine initDirichlet(self, config, object) - use json_module - - class(boundaryEMGeneric), allocatable, intent(inout):: self - type(json_file), intent(inout):: config - character(:), allocatable, intent(in):: object - - end subroutine initDirichlet - - module subroutine initDirichletTime(self, config, object) - use json_module - - class(boundaryEMGeneric), allocatable, intent(inout):: self - type(json_file), intent(inout):: config - character(:), allocatable, intent(in):: object - - end subroutine initDirichletTime - - module function boundaryEMName_to_Index(boundaryName) result(bp) - character(:), allocatable:: boundaryName - integer:: bp - - end function boundaryEMName_to_Index - - end interface - abstract interface ! Apply boundary condition to the load vector for the Poission equation subroutine applyEM_interface(self, vectorF) @@ -1008,6 +998,24 @@ MODULE moduleMesh END TYPE boundaryEMDirichlet + interface + module subroutine initDirichlet(self, config, object) + use json_module + + class(boundaryEMGeneric), allocatable, intent(inout):: self + type(json_file), intent(inout):: config + character(:), allocatable, intent(in):: object + + end subroutine initDirichlet + + module subroutine applyDirichlet(self, vectorF) + class(boundaryEMDirichlet), intent(in):: self + real(8), intent(inout):: vectorF(:) + + end subroutine applyDirichlet + + end interface + TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime real(8):: potential real(8):: timeFactor @@ -1019,11 +1027,14 @@ MODULE moduleMesh END TYPE boundaryEMDirichletTime interface - module subroutine applyDirichlet(self, vectorF) - class(boundaryEMDirichlet), intent(in):: self - real(8), intent(inout):: vectorF(:) + module subroutine initDirichletTime(self, config, object) + use json_module - end subroutine applyDirichlet + class(boundaryEMGeneric), allocatable, intent(inout):: self + type(json_file), intent(inout):: config + character(:), allocatable, intent(in):: object + + end subroutine initDirichletTime module subroutine applyDirichletTime(self, vectorF) class(boundaryEMDirichletTime), intent(in):: self @@ -1047,6 +1058,12 @@ MODULE moduleMesh ! Update the EM boundary models interface + module function boundaryEMName_to_Index(boundaryName) result(bp) + character(:), allocatable:: boundaryName + integer:: bp + + end function boundaryEMName_to_Index + module subroutine boundariesEM_update() end subroutine boundariesEM_update diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index 14fd0ee..8640d24 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -56,6 +56,7 @@ submodule(moduleMesh) boundaryEM use json_module use moduleRefParam, ONLY: Volt_ref, ti_ref use moduleErrors + use moduleOutput, only: path implicit none class(boundaryEMGeneric), allocatable, intent(inout):: self diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index 2e2acfa..04116be 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -342,6 +342,7 @@ submodule(moduleMesh) boundaryParticle ! Print subroutine ionization_print(self, fileID) + use moduleOutput, only: fmtInt implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -473,6 +474,7 @@ submodule(moduleMesh) boundaryParticle ! Print output subroutine quasiNeutrality_print(self, fileID) + use moduleOutput, only: fmtColInt, fmtColReal implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -613,6 +615,7 @@ submodule(moduleMesh) boundaryParticle ! Print subroutine outflowAdaptive_print(self, fileID) + use moduleOutput, only: fmtColInt, fmtColReal implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -691,7 +694,12 @@ submodule(moduleMesh) boundaryParticle ! Writes the output into the Step_XXXXX_boundaryParticles file when the %print procedure is associated module subroutine boundariesParticle_write() use moduleCaseparam, only: timeStep - use moduleOutput, only:fileID_boundaryParticle, formatFileName, informFileCreation + use moduleOutput, only: fileID_boundaryParticle, & + formatFileName, & + informFileCreation,& + boundaryParticleOutput,& + generateFilePath, & + prefix implicit none integer:: b @@ -700,7 +708,7 @@ submodule(moduleMesh) boundaryParticle if (boundaryParticleOutput) then fileName = formatFileName(prefix, 'boundariesParticle', 'csv', timeStep) call informFileCreation(fileName) - open(fileID_boundaryParticle, file = path // folder // '/' // fileName) + open(fileID_boundaryParticle, file = generateFilePath(fileName)) do b = 1, nBoundariesParticle if (associated(boundariesParticle(b)%obj%print)) then diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index d22258a..52e26d3 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -5,14 +5,14 @@ MODULE moduleInject use moduleMesh, only: meshEdgePointer !Generic injection of particles - TYPE:: injectGeneric + TYPE, abstract:: injectGeneric INTEGER:: id CHARACTER(:), ALLOCATABLE:: name REAL(8):: vMod !Velocity (module) 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 + INTEGER:: nParticles !Number of particles to inject each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges type(meshEdgePointer), allocatable:: edges(:) @@ -20,38 +20,125 @@ MODULE moduleInject 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 + procedure(updateInject_interface), pointer, pass:: update => null() + procedure(printInject_interface), pointer, pass:: print => null() CONTAINS - PROCEDURE, PASS:: init => initInject PROCEDURE, PASS:: addParticles END TYPE injectGeneric + abstract interface + ! Update the values of the particle boundary model + subroutine updateInject_interface(self) + import injectGeneric + + class(injectGeneric), intent(inout):: self + + end subroutine updateInject_interface + + ! Write the values of the particle boundary model + subroutine printInject_interface(self, fileID) + import injectGeneric + + class(injectGeneric), intent(inout):: self + integer, intent(in):: fileID + + end subroutine printInject_interface + + end interface + + ! Default type, constant inject of particles + type, extends(injectGeneric):: injectConstant + + end type injectconstant + + type, extends(injectGeneric):: injectQuasiNeutral + real(8):: qSign_incident + + end type injectQuasiNeutral + + ! Wrapper for injects + type:: injectWrapper + class(injectGeneric), allocatable:: obj + + end type injectWrapper + INTEGER:: nInject - TYPE(injectGeneric), ALLOCATABLE:: inject(:) + TYPE(injectWrapper), ALLOCATABLE, target:: injects(:) CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, ps, particlesPerEdge) + SUBROUTINE initInject(self, i, object, config) USE moduleMesh USE moduleRefParam USE moduleConstParam USE moduleSpecies - USE moduleSolver USE moduleErrors + use json_module IMPLICIT NONE - CLASS(injectGeneric), INTENT(inout):: self + CLASS(injectGeneric), allocatable, INTENT(inout):: self INTEGER, INTENT(in):: i - REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) - INTEGER, INTENT(in):: sp, ps, particlesPerEdge + character(:), allocatable, intent(in):: object + type(json_file), intent(inout):: config + character(:), allocatable:: type + logical:: found + REAL(8):: v + real(8), allocatable:: n(:), temperature(:) + INTEGER:: sp, ps, particlesPerEdge + character(:), allocatable:: speciesName integer:: ps_index REAL(8):: tauInject - REAL(8), INTENT(in):: flow - CHARACTER(:), ALLOCATABLE, INTENT(in):: units + REAL(8):: flow + CHARACTER(:), ALLOCATABLE:: units INTEGER:: e REAL(8):: fluxPerStep = 0.D0 - self%id = i + !Type of injection + call config%get(object // '.type', type, found) + if (.not. found) then + ! If no type is found, assume constant injection of particles + type = 'constant' + end if + + ! Assign type to inject and allocates it + select case(type) + case ('constant') + allocate(injectConstant:: self) + + case ('quasiNeutral') + allocate(injectQuasiNeutral:: self) + + self%update => updateQuasiNeutral + self%print => printQuasiNeutral + + case default + call criticalError('No injection type ' // type // ' defined', 'initInject') + end select + + self%id = i + call config%get(object // '.name', self%name, found) + + !Find species + CALL config%get(object // '.species', speciesName, found) + sp = speciesName2Index(speciesName) + + ! Velocity and direction of injection + call config%get(object // '.v', v, found) + call config%get(object // '.n', n, found) + if (.not. found) then + allocate(n(1:3)) + n = 0.D0 + end if + ! Temperature of injection + call config%get(object // '.T', temperature, found) + ! Flow + call config%get(object // '.flow', flow, found) + call config%get(object // '.units', units, found) + CALL config%get(object // '.physicalSurface', ps, found) + particlesPerEdge = 0 + CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) + self%vMod = v / v_ref self%n = n / NORM2(n) self%temperature = temperature / T_ref @@ -124,8 +211,6 @@ MODULE moduleInject END DO - self%nParticles = SUM(self%particlesPerEdge) - ELSE ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight @@ -133,56 +218,235 @@ MODULE moduleInject self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight)) END DO - self%nParticles = SUM(self%particlesPerEdge) - !Rescale weight to match flux self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) END IF + ! Total number of particles to inject + 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') + ! Specific parameters per type + select type(self) + type is(injectQuasiNeutral) + ! Get the sign of the charged species + self%qSign_incident = sign(1.d0, self%species%qm) + + end select + END SUBROUTINE initInject - !Injection of particles - SUBROUTINE doInjects() + ! Update the value of injects + subroutine updateInjects() + implicit none + + integer:: i + + do i = 1, nInject + if (associated(injects(i)%obj%update)) then + call injects(i)%obj%update() + + end if + + end do + + end subroutine updateInjects + + ! Updates the flow in the injection to maintain quasineutrality + subroutine updateQuasiNeutral(self) + use moduleMesh, only: meshEdge, meshCell + implicit none + + class(injectGeneric), intent(inout):: self + integer:: e + class(meshEdge), pointer:: edge + class(meshCell), pointer:: cell + real(8):: Xi(1:3) + real(8):: EF_normal + real(8):: alpha + + select type(self) + type is (injectQuasiNeutral) + do e = 1, self%nEdges + edge => self%edges(e)%obj + + Xi = edge%centerXi() + + if (associated(edge%e1)) then + cell => edge%e1 + + else + cell => edge%e2 + + end if + + ! Projection of EF on the edge normal vector + EF_normal = dot_product(cell%gatherElectricField(Xi), edge%normal) + + ! Correction for this time step + alpha = self%qSign_incident * EF_normal + + ! Adjust the weight of particles to match the new current + self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha + + ! Limit the weight of a macroparticle to 1.0 + self%weightPerEdge(e) = max(self%weightPerEdge(e), 1.d0) + + end do + + end select + + end subroutine updateQuasiNeutral + + !Add particles for the injection + SUBROUTINE addParticles(self) USE moduleSpecies - USE moduleSolver + USE moduleMesh + USE moduleRandom + USE moduleErrors + use moduleList, only: partInj IMPLICIT NONE - INTEGER:: i + CLASS(injectGeneric), INTENT(in):: self + type(particle), pointer:: part + INTEGER:: i, e + INTEGER:: sp + CLASS(meshEdge), POINTER:: edge + + REAL(8):: direction(1:3) + + !Insert particles + !$OMP DO + DO e = 1, self%nEdges + ! Select edge for injection + edge => self%edges(e)%obj + ! Inject particles in edge + DO i = 1, self%particlesPerEdge(e) + allocate(part) + ! Index in the partInj array + !Particle is considered to be inside the domain + part%n_in = .true. + !Random position in edge + part%r = edge%randPos() + !Assign weight to particle. + part%weight = self%weightPerEdge(e) + !Volume associated to the edge: + IF (ASSOCIATED(edge%e1)) THEN + part%cell = edge%e1%n + + ELSEIF (ASSOCIATED(edge%e2)) THEN + part%cell = edge%e2%n + + ELSE + CALL criticalError("No Volume associated to edge", 'addParticles') + + END IF + part%cellColl = edge%eColl%n + sp = self%species%n + + !Assign particle type + part%species => self%species + + if (all(self%n == 0.D0)) then + direction = edge%normal + + else + direction = self%n + + end if + + part%v = 0.D0 + + ! do while(dot_product(partInj(n)%v, direction) <= 0.d0) + part%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) + !If injecting a no-drift distribution and velocity is negative, reflect + if ((self%vMod == 0.D0) .and. & + (dot_product(part%v, direction) <= 0.D0)) then + part%v = - part%v + + end if + + ! end do + + !Obtain natural coordinates of particle in cell + part%Xi = mesh%cells(part%cell)%obj%phy2log(part%r) + + ! Add particle to global list + call partInj%setLock() + call partInj%add(part) + call partInj%unsetLock() + + END DO + + END DO + !$OMP END DO + + END SUBROUTINE addParticles + + ! Writes the value of injects + subroutine writeInjects() + use moduleCaseparam, only: timeStep + use moduleOutput, only: fileID_inject, & + formatFileName, & + informFileCreation,& + injectOutput,& + generateFilePath, & + prefix + implicit none - !$OMP SINGLE - nPartInj = 0 - DO i = 1, nInject - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - nPartInj = nPartInj + inject(i)%nParticles + integer:: i + character(:), allocatable:: fileName - END IF + if (injectOutput) then + fileName = formatFileName(prefix, 'injects', 'csv', timeStep) + call informFileCreation(fileName) + open(fileID_inject, file = generateFilePath(fileName)) - END DO + do i = 1, nInject + if (associated(injects(i)%obj%print)) then + call injects(i)%obj%print(fileID_inject) - IF (ALLOCATED(partInj)) DEALLOCATE(partInj) - ALLOCATE(partInj(1:nPartInj)) - !$OMP END SINGLE + end if - DO i=1, nInject - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - CALL inject(i)%addParticles() + end do - END IF - END DO + close(fileID_inject) - END SUBROUTINE doInjects + end if + end subroutine writeInjects + + ! Write information about quasiNeutral inject + subroutine printQuasiNeutral(self, fileID) + use moduleOutput, only: fmtColInt, fmtColReal + implicit none + + class(injectGeneric), intent(inout):: self + integer, intent(in):: fileID + integer:: e + + write(fileID, '(A)') self%name + write(fileID, '(A,",",A)') '"Edge id"', '"weight of particles"' + do e = 1, self%nEdges + write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%weightPerEdge(e) + + end do + + end subroutine printQuasiNeutral + + ! Inits for different velocity distribution functions SUBROUTINE initVelDistMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -192,7 +456,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -205,101 +469,4 @@ MODULE moduleInject END SUBROUTINE initVelDistDelta - !Add particles for the injection - SUBROUTINE addParticles(self) - USE moduleSpecies - USE moduleSolver - USE moduleMesh - USE moduleRandom - USE moduleErrors - IMPLICIT NONE - - CLASS(injectGeneric), INTENT(in):: self - INTEGER, SAVE:: nMin - INTEGER:: i, e - INTEGER:: n, sp - CLASS(meshEdge), POINTER:: randomEdge - REAL(8):: direction(1:3) - - !Insert particles - !$OMP SINGLE - nMin = 0 - DO i = 1, self%id -1 - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - nMin = nMin + inject(i)%nParticles - - END IF - - END DO - nMin = nMin + 1 - !$OMP END SINGLE - - !$OMP DO - DO e = 1, self%nEdges - ! Select edge for injection - randomEdge => 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 - - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n - - ELSE - CALL criticalError("No Volume associated to edge", 'addParticles') - - END IF - partInj(n)%cellColl = randomEdge%eColl%n - sp = self%species%n - - !Assign particle type - partInj(n)%species => self%species - - if (all(self%n == 0.D0)) then - direction = randomEdge%normal - - else - direction = self%n - - end if - - partInj(n)%v = 0.D0 - - do while(dot_product(partInj(n)%v, direction) <= 0.d0) - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) - !If injecting a no-drift distribution and velocity is negative, reflect - if ((self%vMod == 0.D0) .and. & - (dot_product(partInj(n)%v, direction) <= 0.D0)) then - partInj(n)%v = - partInj(n)%v - - end if - - 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)) - - END DO - - END DO - !$OMP END DO - - END SUBROUTINE addParticles - END MODULE moduleInject diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index b1dafdc..44524c6 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -23,6 +23,7 @@ MODULE moduleList END TYPE listNode + TYPE(listNode):: partInj !Particles comming from injections TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme TYPE(listNode):: partCollisions !Particles created in collisional process TYPE(listNode):: partSurfaces !Particles created in surface interactions diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index ab08f08..8c575da 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -48,11 +48,8 @@ MODULE moduleSpecies !Number of old particles INTEGER:: nPartOld - !Number of injected particles - INTEGER:: nPartInj !Arrays that contain the particles TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles CONTAINS FUNCTION speciesName2Index(speciesName) RESULT(sp) diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index ab2d341..3011396 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -14,6 +14,7 @@ MODULE moduleOutput LOGICAL:: collOutput = .FALSE. LOGICAL:: emOutput = .FALSE. logical:: boundaryParticleOutput = .false. + logical:: injectOutput = .false. ! Prefix for iteration files character(len=*), parameter:: prefix = 'Step' @@ -33,6 +34,7 @@ MODULE moduleOutput integer, parameter:: fileID_output = 20 ! Base id for species/collisions/EM output integer, parameter:: fileID_boundaryParticle = 30 ! Particle boundaries integer, parameter:: fileID_boundaryEM = 31 ! EM boundaries + integer, parameter:: fileID_inject = 32 ! Injects integer, parameter:: fileID_reference = 40 ! Reference values integer, parameter:: fileID_time =50 ! Computation time diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index ed566e9..635c9eb 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -235,23 +235,21 @@ MODULE moduleSolver INTEGER:: nn, n, e INTEGER, SAVE:: nPartNew - INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces + INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) INTEGER:: s !$OMP SECTIONS !$OMP SECTION - nInjIn = 0 - IF (ALLOCATED(partInj)) THEN - nInjIn = COUNT(partInj%n_in) - - END IF - !$OMP SECTION nOldIn = 0 IF (ALLOCATED(partOld)) THEN nOldIn = COUNT(partOld%n_in) END IF + + !$OMP SECTION + nInj = partInj%amount + !$OMP SECTION nWScheme = partWScheme%amount @@ -267,30 +265,14 @@ MODULE moduleSolver !$OMP SINGLE CALL MOVE_ALLOC(partOld, partTemp) - nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces + nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces ALLOCATE(partOld(1:nPartNew)) !$OMP END SINGLE !$OMP SECTIONS - !$OMP SECTION - !Reset particles from injection - nn = 0 - DO n = 1, nPartInj - IF (partInj(n)%n_in) THEN - nn = nn + 1 - partOld(nn) = partInj(n) - IF (nProbes > 0) THEN - CALL doProbes(partOld(nn)) - - END IF - - END IF - - END DO - !$OMP SECTION !Reset particles from previous iteration - nn = nInjIn + nn = 0 DO n = 1, nPartOld IF (partTemp(n)%n_in) THEN nn = nn + 1 @@ -304,19 +286,24 @@ MODULE moduleSolver END DO + !$OMP SECTION + !Reset particles from injection + nn = nOldIn + call resetList(partInj, partOld, nn) + !$OMP SECTION !Reset particles from weighting scheme - nn = nInjIn + nOldIn + nn = nOldIn + nInj CALL resetList(partWScheme, partOld, nn) !$OMP SECTION !Reset particles from collisional process - nn = nInjIn + nOldIn + nWScheme + nn = nOldIn + nInj + nWScheme CALL resetList(partCollisions, partOld, nn) !$OMP SECTION !Reset particles from surface process - nn = nInjIn + nOldIn + nWScheme + nCollisions + nn = nOldIn + nInj + nWScheme + nCollisions CALL resetList(partSurfaces, partOld, nn) !$OMP SECTION @@ -473,6 +460,28 @@ MODULE moduleSolver END SUBROUTINE splitParticle + !Injection of particles + SUBROUTINE doInjects() + USE moduleSpecies + use moduleInject + IMPLICIT NONE + + INTEGER:: i, sp + + DO i=1, nInject + associate(inject => injects(i)%obj) + sp = inject%species%n + IF (solver%pusher(sp)%pushSpecies) THEN + + CALL inject%addParticles() + + END IF + end associate + + END DO + + END SUBROUTINE doInjects + SUBROUTINE updateParticleCell(self, part) USE moduleSpecies USE moduleMesh @@ -520,6 +529,7 @@ MODULE moduleSolver USE moduleSpecies USE moduleCompTime USE moduleProbe + use moduleInject, only: writeInjects USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE @@ -539,6 +549,9 @@ MODULE moduleSolver ! Output of boundaries call boundariesParticle_write() + ! Output of inejcts + call writeInjects() + WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld IF (timeStep == 0) THEN