diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 492eaf5..9494517 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -826,7 +826,6 @@ MODULE moduleInput real(8):: eThreshold !Energy threshold integer:: speciesID, electronSecondaryID character(:), allocatable:: speciesName, crossSection, electronSecondary - integer:: s_incident ! Read models of particles object = 'boundaries.particles' @@ -901,12 +900,7 @@ MODULE moduleInput END IF case('quasiNeutrality') - call config%get(object // '.incident', speciesName, found) - if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary') - - s_incident = speciesName2Index(speciesName) - - call initQuasiNeutrality(self, s_incident) + call initQuasiNeutrality(self) CASE DEFAULT CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary') diff --git a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 index 0cd539b..e243c80 100644 --- a/src/modules/mesh/inout/text/moduleMeshOutputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshOutputText.f90 @@ -137,10 +137,10 @@ module moduleMeshOutputText write(*, "(6X,A15,A)") "Creating file: ", fileName open (fileID, file = path // folder // '/' // fileName) - write(fileID, '(5(A,","),A)') '"Position (m)"', & - '"Density (m^-3)"', & - '"Velocity (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (m s^-1):2"', & - '"Temperature (K)"' + write(fileID, '(5(A,","),A)') 'Position (m)', & + 'Density (m^-3)', & + 'Velocity (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', & + 'Temperature (K)' call writeSpeciesOutput(self, fileID, s) @@ -206,10 +206,10 @@ module moduleMeshOutputText write(*, "(6X,A15,A)") "Creating file: ", fileName open (fileID, file = path // folder // '/' // fileName) - write(fileID, '(8(A,","),A)') '"Position (m)"', & - '"Potential (V)"', & - '"Electric Field (V m^-1):0"', 'Electric Field (V m^-1):1"', 'Electric Field (V m^-1):2"', & - '"Magnetic Field (T):0"', 'Magnetic Field (T):1"', 'Magnetic Field (T):2"' + write(fileID, '(8(A,","),A)') 'Position (m)', & + 'Potential (V)', & + 'Electric Field (V m^-1):0', 'Electric Field (V m^-1):1', 'Electric Field (V m^-1):2', & + 'Magnetic Field (T):0', 'Magnetic Field (T):1', 'Magnetic Field (T):2' call writeEMOutput(self, fileID) @@ -234,23 +234,23 @@ module moduleMeshOutputText fileIDDeviation = 67 do s = 1, nSpecies - fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv') + fileNameMean = formatFileName('Average_mean', species(s)%obj%name, 'csv', timeStep) write(*, "(6X,A15,A)") "Creating file: ", fileNameMean open (fileIDMean, file = path // folder // '/' // fileNameMean) - 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"', & - '"Temperature, mean (K)"' + 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', & + 'Temperature, mean (K)' - fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv') + fileNameDeviation = formatFileName('Average_deviation', species(s)%obj%name, 'csv', timeStep) write(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation open (fileIDDeviation, file = path // folder // '/' // fileNameDeviation) - 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"', & - '"Temperature, deviation (K)"' + 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', & + 'Temperature, deviation (K)' call writeAverage(self, fileIDMean, fileIDDeviation, s) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 25d7175..9133cca 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -675,9 +675,8 @@ MODULE moduleMesh end subroutine initIonization - module subroutine initQuasiNeutrality(boundary, s_incident) + module subroutine initQuasiNeutrality(boundary) class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - integer, intent(in):: s_incident end subroutine initQuasiNeutrality @@ -763,7 +762,7 @@ MODULE moduleMesh ! Ensures quasi-neutrality by changing the reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality - real(8), allocatable:: alpha(:) ! Reflection parameter + real(8):: alpha ! Reflection parameter integer:: s_incident ! species index of the incident species type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index c51c26d..040f633 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -95,22 +95,17 @@ submodule(moduleMesh) boundaryParticle END SUBROUTINE initIonization - module subroutine initQuasiNeutrality(boundary, s_incident) + module subroutine initQuasiNeutrality(boundary) implicit none class(boundaryParticleGeneric), allocatable, intent(inout):: boundary - integer, intent(in):: s_incident allocate(boundaryQuasiNeutrality:: boundary) select type(boundary) type is(boundaryQuasiNeutrality) - allocate(boundary%edges(0)) - - boundary%s_incident = s_incident - - allocate(boundary%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here boundary%alpha = 0.d0 + allocate(boundary%edges(0)) boundary%update => quasiNeutrality_update @@ -338,7 +333,7 @@ submodule(moduleMesh) boundaryParticle class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - if (random() <= self%alpha(edge%n)) then + if (random() <= self%alpha) then call genericReflection(edge, part) else @@ -352,54 +347,12 @@ submodule(moduleMesh) boundaryParticle implicit none class(boundaryParticleGeneric), intent(inout):: self - integer:: e, n, s - integer, allocatable:: nodes(:) - class(meshEdge), pointer:: edge - real(8), allocatable:: density_nodes(:) - class(meshNode), pointer:: node - real(8):: density_incident, density_rest - real(8):: alpha select type(self) type is(boundaryQuasiNeutrality) - do e = 1, size(self%edges) - edge => mesh%edges(e)%obj + self%alpha = 0.1d0 - 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(n)%obj - density_nodes(n) = node%output(s)%den - - end do - - if (s == self%s_incident) then - density_incident = edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes) - - else - density_rest = density_rest + edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes) - - end if - - end do - - ! Correction for this time step - alpha = 1.d0 - density_incident/density_rest - - ! Apply correction with a factor of 0.1 to avoid fast changes - self%alpha(edge%n) = self%alpha(edge%n) + 0.1d0 * alpha - - ! Limit alpha between 0 and 1 - self%alpha(edge%n) = min(self%alpha(edge%n), 1.d0) - self%alpha(edge%n) = max(self%alpha(edge%n), 0.d0) - - deallocate(density_nodes) - - end do + print *, self%alpha end select