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)
+            call initQuasiNeutrality(self)
-            if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary')
-
-            s_incident = speciesName2Index(speciesName)
-
-            call initQuasiNeutrality(self, s_incident)
 
           CASE DEFAULT
             CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')

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)"', &
+        write(fileID, '(5(A,","),A)') 'Position (m)', &
-                                      '"Density (m^-3)"', &
+                                      'Density (m^-3)', &
-                                      '"Velocity (m s^-1):0"', 'Velocity (m s^-1):1"', 'Velocity (m s^-1):2"', &
+                                      'Velocity (m s^-1):0', 'Velocity (m s^-1):1', 'Velocity (m s^-1):2', &
-                                      '"Temperature (K)"'
+                                      '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)"', &
+        write(fileID, '(8(A,","),A)') 'Position (m)', &
-                                      '"Potential (V)"', &
+                                      'Potential (V)', &
-                                      '"Electric Field (V m^-1):0"', 'Electric Field (V m^-1):1"', 'Electric Field (V m^-1):2"', &
+                                      '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"'
+                                      '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)"', &
+        write(fileIDMean, '(5(A,","),A)') 'Position (m)', &
-                                          '"Density, mean (m^-3)"', &
+                                          '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 (m s^-1):1', 'Velocity (m s^-1):2', &
-                                          '"Temperature, mean (K)"'
+                                          '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)"', &
+        write(fileIDDeviation, '(5(A,","),A)') 'Position (m)', &
-                                               '"Density, deviation (m^-3)"', &
+                                               '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 (m s^-1):1', 'Velocity (m s^-1):2', &
-                                               '"Temperature, deviation (K)"'
+                                               'Temperature, deviation (K)'
 
         call writeAverage(self, fileIDMean, fileIDDeviation, s)
 

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

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)
+        self%alpha = 0.1d0
-          edge => mesh%edges(e)%obj
 
-          density_incident = 0.d0
+        print *, self%alpha
-          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
 
       end select