diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 48c3657..feee89a 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -72,7 +72,7 @@ PROGRAM fpakc call boundariesEM_update() ! Update injects - call updateInjects + 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..f4dcc46 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) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 4f4abae..bcf4dcc 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,50 +1428,17 @@ 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 - character(:), allocatable:: type - INTEGER:: physicalSurface - INTEGER:: particlesPerEdge - INTEGER:: sp CALL config%info('inject', found, n_children = nInject) - ALLOCATE(inject(1:nInject)) + ALLOCATE(injects(1:nInject)) nPartInj = 0 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 // '.type', type, found) - if (.not. found) then - ! If no type is found, assume constant injection of particles - type = 'constant' - end if - 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, type, sp, physicalSurface, particlesPerEdge) - - CALL readVelDistr(config, inject(i), object) + CALL readVelDistr(config, injects(i)%obj, object) END DO @@ -1528,7 +1497,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/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..168f519 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 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 427b3ce..39de1ad 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -5,7 +5,7 @@ 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) @@ -14,46 +14,131 @@ MODULE moduleInject LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n INTEGER:: nParticles !Number of particles to introduce each time step CLASS(speciesGeneric), POINTER:: species !Species of injection - character(:), allocatable:: type INTEGER:: nEdges type(meshEdgePointer), allocatable:: edges(:) INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge REAL(8):: surface ! Total surface of injection TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction + 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 + + end type injectQuasiNeutral + + ! Wrapper for injects + type:: injectWrapper + class(injectGeneric), allocatable:: obj + + end type injectWrapper + INTEGER:: nInject - TYPE(injectGeneric), ALLOCATABLE, target:: inject(:) + TYPE(injectWrapper), ALLOCATABLE, target:: injects(:) CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, temperature, flow, units, type, 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 - character(:), allocatable, intent(in):: type + 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 @@ -145,15 +230,89 @@ MODULE moduleInject !Scale particles for different species steps IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') - ! Assign type to inject - select case(type) - case ('constant', 'quasiNeutral') - self%type = type - case default - call criticalError('No injection type ' // type // ' defined', 'initInject') + 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, meshNode, mesh, qSpecies + implicit none + + class(injectGeneric), intent(inout):: self + integer:: e, s, n + integer, allocatable:: nodes(:) + class(meshEdge), pointer:: edge + real(8), allocatable:: density_nodes(:) + class(meshNode), pointer:: node + real(8):: den_center + real(8):: density_incident, density_rest + real(8):: alpha + + select type(self) + type is (injectQuasiNeutral) + do e = 1, self%nEdges + edge => self%edges(e)%obj + + density_incident = 0.d0 + density_rest = 0.d0 + + nodes = edge%getNodes(edge%nNodes) + allocate(density_nodes(1:edge%nNodes)) + do s = 1, nSpecies + do n = 1, edge%nNodes + node => mesh%nodes(nodes(n))%obj + + density_nodes(n) = node%output(s)%den + + end do + + ! Gather the density at the edge center and multiply by the species charge + den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes) + + if (s == self%species%n) then + density_incident = den_center + + else + density_rest = density_rest + den_center + + end if + + end do + + ! Correction for this time step + if (density_rest > 1.0d-10) then + ! If there is a rest population, correct + alpha = 1.d0 + density_incident/density_rest + + else + ! If not, alpha is assumed unchaged + alpha = 0.d0 + + end if + + ! Adjust the weight of particles to match the new current + self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha + + end do + end select - END SUBROUTINE initInject + end subroutine updateQuasiNeutral !Injection of particles SUBROUTINE doInjects() @@ -166,8 +325,8 @@ MODULE moduleInject !$OMP SINGLE nPartInj = 0 DO i = 1, nInject - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - nPartInj = nPartInj + inject(i)%nParticles + IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN + nPartInj = nPartInj + injects(i)%obj%nParticles END IF @@ -178,112 +337,14 @@ MODULE moduleInject !$OMP END SINGLE DO i=1, nInject - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - CALL inject(i)%addParticles() + IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN + CALL injects(i)%obj%addParticles() END IF END DO END SUBROUTINE doInjects - ! Update the value of injects - subroutine updateInjects() - use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies - implicit none - - integer:: i, e, s, n - class(injectGeneric), pointer:: self - integer, allocatable:: nodes(:) - class(meshEdge), pointer:: edge - real(8), allocatable:: density_nodes(:) - class(meshNode), pointer:: node - real(8):: den_center - real(8):: density_incident, density_rest - real(8):: alpha - - do i = 1, nInject - self => inject(i) - select case(inject(i)%type) - case ('quasiNeutral') - do e = 1, self%nEdges - edge => self%edges(e)%obj - - density_incident = 0.d0 - density_rest = 0.d0 - - nodes = edge%getNodes(edge%nNodes) - allocate(density_nodes(1:edge%nNodes)) - do s = 1, nSpecies - do n = 1, edge%nNodes - node => mesh%nodes(nodes(n))%obj - - density_nodes(n) = node%output(s)%den - - end do - - ! Gather the density at the edge center and multiply by the species charge - den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes) - - if (s == self%species%n) then - density_incident = den_center - - else - density_rest = density_rest + den_center - - end if - - end do - - ! Correction for this time step - if (density_rest > 1.0d-10) then - ! If there is a rest population, correct - alpha = 1.d0 + density_incident/density_rest - - else - ! If not, alpha is assumed unchaged - alpha = 0.d0 - - end if - - ! Adjust the weight of particles to match the new current - self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-2 * alpha - - end do - - end select - end do - - end subroutine updateInjects - - 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)) - - END SUBROUTINE initVelDistMaxwellian - - SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m) - IMPLICIT NONE - - CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: temperature, m - - velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) - - END SUBROUTINE initVelDistHalfMaxwellian - - SUBROUTINE initVelDistDelta(velDist) - IMPLICIT NONE - - CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - - velDist = velDistDelta() - - END SUBROUTINE initVelDistDelta - !Add particles for the injection SUBROUTINE addParticles(self) USE moduleSpecies @@ -304,8 +365,8 @@ MODULE moduleInject !$OMP SINGLE nMin = 0 DO i = 1, self%id -1 - IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN - nMin = nMin + inject(i)%nParticles + IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN + nMin = nMin + injects(i)%obj%nParticles END IF @@ -381,4 +442,85 @@ MODULE moduleInject 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 + + integer:: i + character(:), allocatable:: fileName + + if (injectOutput) then + fileName = formatFileName(prefix, 'injects', 'csv', timeStep) + call informFileCreation(fileName) + open(fileID_inject, file = generateFilePath(fileName)) + + do i = 1, nInject + if (associated(injects(i)%obj%print)) then + call injects(i)%obj%print(fileID_inject) + + 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(2.d0*temperature/m)) + + END SUBROUTINE initVelDistMaxwellian + + SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m) + IMPLICIT NONE + + CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist + REAL(8), INTENT(in):: temperature, m + + velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + + END SUBROUTINE initVelDistHalfMaxwellian + + SUBROUTINE initVelDistDelta(velDist) + IMPLICIT NONE + + CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist + + velDist = velDistDelta() + + END SUBROUTINE initVelDistDelta + END MODULE moduleInject 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..f951839 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -520,6 +520,7 @@ MODULE moduleSolver USE moduleSpecies USE moduleCompTime USE moduleProbe + use moduleInject, only: writeInjects USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE @@ -539,6 +540,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