diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index f8ab920..c3b66ab 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 75f0106..1b4e258 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -401,9 +401,6 @@ 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. @@ -691,14 +688,6 @@ 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 cb68754..921d4ac 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -1,39 +1,15 @@ ! FPAKC main program PROGRAM fpakc - 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 + USE moduleCompTime + USE moduleCaseParam + USE moduleInput + USE moduleInject + USE moduleSolver + USE moduleMesh + USE moduleProbe + USE moduleErrors + USE OMP_LIB + IMPLICIT NONE ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 @@ -95,9 +71,6 @@ 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 a1bc3fc..1b8ec67 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 + real(8):: v1, v2, Rsquare v1 = 0.d0 do while (v1 <= 0.d0) @@ -67,20 +67,19 @@ 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):: v1 + REAL(8):: rnd + REAL(8):: x rnd = 0.D0 - v1 = 0.D0 - do while (v1 == 0.D0) - v1 = random() - - end do + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + END DO - rnd = sqrt(-2.d0*log(v1)) + rnd = DSQRT(-DLOG(x)) END FUNCTION randomHalfMaxwellian diff --git a/src/modules/common/velocityDistribution.f90 b/src/modules/common/velocityDistribution.f90 index 359eebd..5b6c342 100644 --- a/src/modules/common/velocityDistribution.f90 +++ b/src/modules/common/velocityDistribution.f90 @@ -53,11 +53,12 @@ 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() + v = self%vTh*randomMaxwellian()/sqrt(2.d0) end function randomVelMaxwellian diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index eb8bde5..8b7aaff 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -493,7 +493,6 @@ 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 @@ -629,7 +628,6 @@ MODULE moduleInput USE moduleMesh USE moduleCaseParam USE moduleRefParam - use moduleOutput, only: collOutput USE OMP_LIB USE json_module IMPLICIT NONE @@ -1428,19 +1426,53 @@ 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(injects(1:nInject)) + ALLOCATE(inject(1:nInject)) + nPartInj = 0 DO i = 1, nInject WRITE(iString, '(i2)') i object = 'inject(' // trim(iString) // ')' - CALL initInject(injects(i)%obj, i, object, config) + !Find species + CALL config%get(object // '.species', speciesName, found) + sp = speciesName2Index(speciesName) - CALL readVelDistr(config, injects(i)%obj, object) + 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) END DO + !Allocate array for injected particles + IF (nPartInj > 0) THEN + ALLOCATE(partInj(1:nPartInj)) + + END IF + END SUBROUTINE readInject SUBROUTINE readAverage(config) @@ -1490,7 +1522,7 @@ MODULE moduleInput IMPLICIT NONE TYPE(json_file), INTENT(inout):: config - class(injectGeneric), INTENT(inout):: inj + TYPE(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 5789276..9f62c85 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: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o +solver.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 acd4189..0ca61ad 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -1,5 +1,4 @@ MODULE moduleMeshOutputGmsh2 - use moduleOutput, only: prefix, informFileCreation, formatFileName, generateFilePath CONTAINS !Header for mesh format @@ -17,7 +16,7 @@ MODULE moduleMeshOutputGmsh2 !Node data subroutines !Header SUBROUTINE writeGmsh2HeaderNodeData(fileID, title, iteration, time, dimensions, nNodes) - use moduleOutput, only: fmtInt, fmtReal + use moduleOutput IMPLICIT NONE INTEGER, INTENT(in):: fileID @@ -51,7 +50,7 @@ MODULE moduleMeshOutputGmsh2 !Element data subroutines !Header SUBROUTINE writeGmsh2HeaderElementData(fileID, title, iteration, time, dimensions, nVols) - use moduleOutput, only: fmtInt, fmtReal + use moduleOutput IMPLICIT NONE INTEGER, INTENT(in):: fileID @@ -87,7 +86,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleMesh USE moduleRefParam USE moduleSpecies - use moduleOutput, only: outputFormat, calculateOutput, fmtReal + USE moduleOutput USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE @@ -142,7 +141,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleRefParam USE moduleCaseParam USE moduleCollisions - use moduleOutput, only: collOutput, fmtInt + USE moduleOutput IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self @@ -200,7 +199,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleMesh USE moduleRefParam USE moduleCaseParam - use moduleOutput, only: emOutput + USE moduleOutput IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -249,6 +248,7 @@ 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 74ac0ea..7f95788 100644 --- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 @@ -1,14 +1,9 @@ module moduleMeshOutputText - use moduleOutput, only: prefix, & - informFileCreation, & - formatFileName, & - generateFilePath - contains subroutine writeSpeciesOutput(self, fileID, speciesIndex) use moduleMesh - use moduleOutput, only: calculateOutput, outputFormat + use moduleOutput use moduleRefParam, only: L_ref implicit none @@ -93,6 +88,7 @@ module moduleMeshOutputText speciesIndex) use moduleMesh + use moduleOutput use moduleAverage use moduleRefParam, only: L_ref implicit none @@ -155,8 +151,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 @@ -194,7 +190,6 @@ module moduleMeshOutputText subroutine printEMText(self) use moduleMesh use moduleCaseParam, only: timeStep - use moduleOutput, only: emOutput implicit none class(meshParticles), intent(in):: self @@ -241,7 +236,7 @@ module moduleMeshOutputText write(fileIDMean, '(5(A,","),A)') '"Position (m)"', & '"Density, mean (m^-3)"', & - '"Velocity, mean (m s^-1):0"', '"Velocity, mean (m s^-1):1"', '"Velocity, mean (m s^-1):2"', & + '"Velocity, mean (m s^-1):0"', '"Velocity (m s^-1):1"', '"Velocity (m s^-1):2"', & '"Temperature, mean (K)"' fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv') @@ -250,7 +245,7 @@ module moduleMeshOutputText write(fileIDDeviation, '(5(A,","),A)') '"Position (m)"', & '"Density, deviation (m^-3)"', & - '"Velocity, deviation (m s^-1):0"', 'Velocity, deviation (m s^-1):1"', 'Velocity, deviation (m s^-1):2"', & + '"Velocity, deviation (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (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 8058896..3462940 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -1,11 +1,9 @@ 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 @@ -79,7 +77,6 @@ MODULE moduleMeshOutputVTU SUBROUTINE writeMesh(self, fileID) USE moduleMesh USE moduleRefParam - use moduleOutput, only: fmtReal, fmtInt IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self @@ -162,7 +159,6 @@ MODULE moduleMeshOutputVTU SUBROUTINE writeCollOutput(self,fileID) USE moduleMesh USE moduleCollisions - use moduleOutput, only: fmtInt IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self @@ -189,7 +185,6 @@ MODULE moduleMeshOutputVTU SUBROUTINE writeEM(self, fileID) USE moduleMesh USE moduleRefParam - use moduleOutput, only: fmtReal IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -386,7 +381,6 @@ MODULE moduleMeshOutputVTU SUBROUTINE printEMVTU(self) USE moduleMesh USE moduleCaseParam, ONLY: timeStep - use moduleOutput, only: emOutput IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -420,6 +414,7 @@ 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 6cc48ee..7ce7565 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, only: outputNode, emNode + USE moduleOutput USE moduleCollisions use moduleSpecies, only: nSpecies @@ -663,6 +663,50 @@ 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) @@ -701,18 +745,6 @@ 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 @@ -720,18 +752,6 @@ 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 @@ -739,18 +759,6 @@ 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 @@ -758,18 +766,6 @@ 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) @@ -779,24 +775,6 @@ 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. @@ -813,30 +791,6 @@ 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 @@ -847,24 +801,6 @@ 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(:) @@ -876,11 +812,68 @@ MODULE moduleMesh end type boundaryOutflowAdaptive interface - module subroutine outflowAdaptive_init(boundary, s_incident) - class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - integer, intent(in):: s_incident + module subroutine reflection_apply(self, edge, part) + use moduleSpecies - end subroutine outflowAdaptive_init + 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 module subroutine outflowAdaptive_apply(self, edge, part) use moduleSpecies @@ -891,10 +884,7 @@ MODULE moduleMesh end subroutine outflowAdaptive_apply - end interface - - ! Generic basic boundary conditions to use internally in the code - interface + ! Generic basic boundary conditions to use internally in the code module subroutine genericReflection(edge, part) use moduleSpecies @@ -937,13 +927,6 @@ 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() @@ -969,6 +952,33 @@ 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) @@ -998,24 +1008,6 @@ 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 @@ -1027,14 +1019,11 @@ MODULE moduleMesh END TYPE boundaryEMDirichletTime interface - module subroutine initDirichletTime(self, config, object) - use json_module + module subroutine applyDirichlet(self, vectorF) + class(boundaryEMDirichlet), intent(in):: self + real(8), intent(inout):: vectorF(:) - class(boundaryEMGeneric), allocatable, intent(inout):: self - type(json_file), intent(inout):: config - character(:), allocatable, intent(in):: object - - end subroutine initDirichletTime + end subroutine applyDirichlet module subroutine applyDirichletTime(self, vectorF) class(boundaryEMDirichletTime), intent(in):: self @@ -1058,12 +1047,6 @@ 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 8640d24..14fd0ee 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -56,7 +56,6 @@ 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 04116be..2e2acfa 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -342,7 +342,6 @@ submodule(moduleMesh) boundaryParticle ! Print subroutine ionization_print(self, fileID) - use moduleOutput, only: fmtInt implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -474,7 +473,6 @@ submodule(moduleMesh) boundaryParticle ! Print output subroutine quasiNeutrality_print(self, fileID) - use moduleOutput, only: fmtColInt, fmtColReal implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -615,7 +613,6 @@ submodule(moduleMesh) boundaryParticle ! Print subroutine outflowAdaptive_print(self, fileID) - use moduleOutput, only: fmtColInt, fmtColReal implicit none class(boundaryParticleGeneric), intent(inout):: self @@ -694,12 +691,7 @@ 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,& - boundaryParticleOutput,& - generateFilePath, & - prefix + use moduleOutput, only:fileID_boundaryParticle, formatFileName, informFileCreation implicit none integer:: b @@ -708,7 +700,7 @@ submodule(moduleMesh) boundaryParticle if (boundaryParticleOutput) then fileName = formatFileName(prefix, 'boundariesParticle', 'csv', timeStep) call informFileCreation(fileName) - open(fileID_boundaryParticle, file = generateFilePath(fileName)) + open(fileID_boundaryParticle, file = path // folder // '/' // 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 52e26d3..d22258a 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, abstract:: injectGeneric + TYPE:: 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 inject each time step + INTEGER:: nParticles !Number of particles to introduce each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges type(meshEdgePointer), allocatable:: edges(:) @@ -20,125 +20,38 @@ 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(injectWrapper), ALLOCATABLE, target:: injects(:) + TYPE(injectGeneric), ALLOCATABLE:: inject(:) CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, object, config) + SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, ps, particlesPerEdge) USE moduleMesh USE moduleRefParam USE moduleConstParam USE moduleSpecies + USE moduleSolver USE moduleErrors - use json_module IMPLICIT NONE - CLASS(injectGeneric), allocatable, INTENT(inout):: self + CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i - 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 + REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) + INTEGER, INTENT(in):: sp, ps, particlesPerEdge integer:: ps_index REAL(8):: tauInject - REAL(8):: flow - CHARACTER(:), ALLOCATABLE:: units + REAL(8), INTENT(in):: flow + CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e REAL(8):: fluxPerStep = 0.D0 - !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%id = i self%vMod = v / v_ref self%n = n / NORM2(n) self%temperature = temperature / T_ref @@ -211,6 +124,8 @@ MODULE moduleInject END DO + self%nParticles = SUM(self%particlesPerEdge) + ELSE ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight @@ -218,235 +133,56 @@ 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 - ! 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) + !Injection of particles + SUBROUTINE doInjects() USE moduleSpecies - USE moduleMesh - USE moduleRandom - USE moduleErrors - use moduleList, only: partInj + USE moduleSolver IMPLICIT NONE - CLASS(injectGeneric), INTENT(in):: self - type(particle), pointer:: part - INTEGER:: i, e - INTEGER:: sp - CLASS(meshEdge), POINTER:: edge + INTEGER:: i + + !$OMP SINGLE + nPartInj = 0 + DO i = 1, nInject + IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN + nPartInj = nPartInj + inject(i)%nParticles - 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 IF END DO - !$OMP END DO - END SUBROUTINE addParticles + IF (ALLOCATED(partInj)) DEALLOCATE(partInj) + ALLOCATE(partInj(1:nPartInj)) + !$OMP END SINGLE - ! Writes the value of injects - subroutine writeInjects() - use moduleCaseparam, only: timeStep - use moduleOutput, only: fileID_inject, & - formatFileName, & - informFileCreation,& - injectOutput,& - generateFilePath, & - prefix - implicit none - - integer:: i - character(:), allocatable:: fileName + DO i=1, nInject + IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN + CALL inject(i)%addParticles() - if (injectOutput) then - fileName = formatFileName(prefix, 'injects', 'csv', timeStep) - call informFileCreation(fileName) - open(fileID_inject, file = generateFilePath(fileName)) + END IF + END DO - do i = 1, nInject - if (associated(injects(i)%obj%print)) then - call injects(i)%obj%print(fileID_inject) + END SUBROUTINE doInjects - end if - - end do - - close(fileID_inject) - - 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(temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -456,7 +192,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -469,4 +205,101 @@ 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 44524c6..b1dafdc 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -23,7 +23,6 @@ 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 8c575da..ab08f08 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -48,8 +48,11 @@ 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 3011396..ab2d341 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -14,7 +14,6 @@ MODULE moduleOutput LOGICAL:: collOutput = .FALSE. LOGICAL:: emOutput = .FALSE. logical:: boundaryParticleOutput = .false. - logical:: injectOutput = .false. ! Prefix for iteration files character(len=*), parameter:: prefix = 'Step' @@ -34,7 +33,6 @@ 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 635c9eb..ed566e9 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -235,21 +235,23 @@ MODULE moduleSolver INTEGER:: nn, n, e INTEGER, SAVE:: nPartNew - INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces + INTEGER, SAVE:: nInjIn, nOldIn, 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 @@ -265,14 +267,30 @@ MODULE moduleSolver !$OMP SINGLE CALL MOVE_ALLOC(partOld, partTemp) - nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces + nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces ALLOCATE(partOld(1:nPartNew)) !$OMP END SINGLE !$OMP SECTIONS !$OMP SECTION - !Reset particles from previous iteration + !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 DO n = 1, nPartOld IF (partTemp(n)%n_in) THEN nn = nn + 1 @@ -286,24 +304,19 @@ 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 = nOldIn + nInj + nn = nInjIn + nOldIn CALL resetList(partWScheme, partOld, nn) !$OMP SECTION !Reset particles from collisional process - nn = nOldIn + nInj + nWScheme + nn = nInjIn + nOldIn + nWScheme CALL resetList(partCollisions, partOld, nn) !$OMP SECTION !Reset particles from surface process - nn = nOldIn + nInj + nWScheme + nCollisions + nn = nInjIn + nOldIn + nWScheme + nCollisions CALL resetList(partSurfaces, partOld, nn) !$OMP SECTION @@ -460,28 +473,6 @@ 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 @@ -529,7 +520,6 @@ MODULE moduleSolver USE moduleSpecies USE moduleCompTime USE moduleProbe - use moduleInject, only: writeInjects USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE @@ -549,9 +539,6 @@ 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