Merge branch 'feature/average' into 'development'

Average scheme

See merge request JorgeGonz/fpakc!29

doc/user-manual/bibliography.bib

@@ -51,4 +51,15 @@
   howpublished = {\url{https://gmsh.info/}},
 }
 
+@Article{welford1962note,
+  author    = {Welford, BP},
+  journal   = {Technometrics},
+  title     = {Note on a method for calculating corrected sums of squares and products},
+  year      = {1962},
+  number    = {3},
+  pages     = {419--420},
+  volume    = {4},
+  publisher = {Taylor \& Francis},
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}

doc/user-manual/fpakc_UserManual.tex

@@ -1,4 +1,4 @@
-\documentclass[10pt,a4paper,oneside]{book}
+\documentclass[10pt,a4paper,twoside]{book}
 \usepackage[latin1]{inputenc}
 \usepackage{amsmath}
 \usepackage{amsfonts}
@@ -22,6 +22,7 @@
     Author   = {Jorge Gonzalez}
   }
 }
+\usepackage[margin=1in]{geometry} % Reduces margins of the document
 % Allows breaking of URL in bibliography.
 \setcounter{biburllcpenalty}{7000}
 \setcounter{biburlucpenalty}{8000}
@@ -260,6 +261,18 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Electromagnetic field}
       WIP.
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    \section{Average scheme}
+      Particle-in-cell codes has an intrinsic statistical noise associated with them.
+      Although this can be reduced by increasing the number of particles, this also increases the CPU requirements of the case.
+      
+      It is quite common that most cases reach a quasi-steady state after a number of iterations and time-average results can be obtained after to improve analysis, plotting and restarting the case using these time-average results as new species backgrounds.
+      Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the amount of iterations printed.
+      
+      \Gls{fpakc} implements a simple average scheme that, after a start time provided by the user, scores a mean and standard deviation of all the main species properties, and the electromagnetic field.
+      This scheme is based on the Welford's online algorithm~\cite{welford1962note}.
+      The averaged data is written in the same format as the input mesh at the end of the simulation.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
   \chapter{Installation}
@@ -642,7 +655,7 @@ make
           \item \textbf{initialTime}: Real.
                                       Units of $\unit{s}$.
                                       Initial simulation time.
-                                      If no value is provided, the initial time is set to $\unit[0]{s}$.
+                                      If no value is provided, the initial time is set to $\unit[0.0]{s}$.
           \item \textbf{pusher}: Character.
                                  Array dimension 'number of species'.
                                  Indicates the type of pusher used for each species:
@@ -684,6 +697,17 @@ make
                                                             File must be located at \textbf{output.path}.
                                   \end{itemize}
         \end{itemize}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      \subsection{average}
+        This object determines the use of an average scheme.
+        If this object exists in the input file, average will be written at the end of the simulation.
+        Acceptable values are:
+        \begin{itemize}
+          \item \textbf{startTime}: Real.
+                                    Units in $\unit{s}$.
+                                    Simulation physical time in which average scheme will start to compute the mean and standard validation.
+                                    If no value is provided, the initial time is set to $\unit[0.0]{s}$.
+        \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       \subsection{interactions}\label{ssec:input_interactions}
         This object determine the different interactions among species.

src/fpakc.f90

@@ -114,6 +114,8 @@ PROGRAM fpakc
     !$OMP SINGLE
     tEMField = omp_get_wtime() - tEMField
 
+    CALL doAverage(t)
+
     tStep = omp_get_wtime() - tStep
 
     !Output data

src/makefile

@@ -4,7 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
 			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
-					$(OBJDIR)/moduleProbe.o \
+					$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \
 			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
 			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \

src/modules/makefile

@@ -7,7 +7,7 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
 
 all: $(OBJS)
 
-mesh.o: moduleCollisions.o moduleBoundary.o
+mesh.o: moduleCollisions.o moduleBoundary.o moduleAverage.o
 	$(MAKE) -C mesh all
 
 moduleCollisions.o: moduleList.o moduleMath.o moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
@@ -25,7 +25,7 @@ moduleList.o: moduleConstParam.o moduleErrors.o moduleCaseParam.o moduleSpecies.
 moduleOutput.o: moduleMath.o moduleRefParam.o moduleOutput.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleSolver.o: moduleProbe.o moduleEM.o moduleSolver.f90
+moduleSolver.o: moduleProbe.o moduleEM.o moduleAverage.o moduleSolver.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleProbe.o: mesh.o moduleProbe.f90

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -260,10 +260,17 @@ MODULE moduleMesh1DRad
       Xii = 0.D0
       fPsi = self%fPsi(Xii)
       detJ = self%detJac(Xii)
+      !Computes total volume of the cell
       r    = DOT_PRODUCT(fPsi, self%r)
       l = 2.D0*detJ
-      self%volume  =      r*l
+      self%volume  = r*l
-      self%arNodes = fPsi*r*l
+      !Computes volume per node
+      xi = (/-5.D-1,  0.D0, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      self%arNodes(1) = fPsi(1)*r*l
+      xi = (/ 5.D-1,  0.D0, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      self%arNodes(2) = fPsi(2)*r*l
 
     END SUBROUTINE areaRad
 

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -249,7 +249,7 @@ MODULE moduleMesh2DCyl
       dr  = self%r(2) - self%r(1)
       dz  = self%z(2) - self%z(1)
       IF (dr /= 0.D0) THEN
-        r(2) = dr*DSQRT(rnd) + self%r(1)
+        r(2) = dr * DSQRT(rnd) + self%r(1)
         r(1) = dz * (r(2) - self%r(1))/dr + self%z(1)
 
       ELSE
@@ -318,9 +318,22 @@ MODULE moduleMesh2DCyl
       xi = 0.D0
       detJ = self%detJac(xi)*PI8 !4*2*pi
       fPsi = self%fPsi(xi)
+      !Computes total volume of the cell
       r = DOT_PRODUCT(fPsi,self%r)
-      self%volume  =      r*detJ
+      self%volume  = r*detJ
-      self%arNodes = fPsi*r*detJ
+      !Computes volume per node
+      xi = (/-5.D-1, -5.D-1, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      self%arNodes(1) = fPsi(1)*r*detJ
+      xi = (/ 5.D-1, -5.D-1, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      self%arNodes(2) = fPsi(2)*r*detJ
+      xi = (/ 5.D-1,  5.D-1, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      self%arNodes(3) = fPsi(3)*r*detJ
+      xi = (/-5.D-1,  5.D-1, 0.D0/)
+      r = DOT_PRODUCT(self%fPsi(xi),self%r)
+      self%arNodes(4) = fPsi(4)*r*detJ
 
     END SUBROUTINE areaQuad
 
@@ -701,8 +714,10 @@ MODULE moduleMesh2DCyl
       xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
       detJ  = self%detJac(xi)*PI !2PI*1/2
       fPsi  = self%fPsi(xi)
+      !Computes total volume of the cell
       r = DOT_PRODUCT(fPsi,self%r)
       self%volume  =      r*detJ
+      !Computes volume per node
       self%arNodes = fPsi*r*detJ
 
     END SUBROUTINE areaTria

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -13,9 +13,10 @@ MODULE moduleMeshInputGmsh2
       IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollGmsh2
       SELECT TYPE(self)
       TYPE IS(meshParticles)
-        self%printOutput => printOutputGmsh2
+        self%printOutput  => printOutputGmsh2
-        self%printEM     => printEMGmsh2
+        self%printEM      => printEMGmsh2
-        self%readInitial => readInitialGmsh2
+        self%readInitial  => readInitialGmsh2
+        self%printAverage => printAverageGmsh2
 
       END SELECT
       self%readMesh    => readGmsh2

src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90

@@ -1,6 +1,84 @@
 MODULE moduleMeshOutputGmsh2
 
   CONTAINS
+    !Header for mesh format
+    SUBROUTINE writeGmsh2HeaderMesh(fileID)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: fileID
+
+      WRITE(fileID, "(A)") '$MeshFormat'
+      WRITE(fileID, "(A)") '2.2 0 8'
+      WRITE(fileID, "(A)") '$EndMeshFormat'
+
+    END SUBROUTINE writeGmsh2HeaderMesh
+
+    !Node data subroutines
+    !Header
+    SUBROUTINE writeGmsh2HeaderNodeData(fileID, title, iteration, time, dimensions, nNodes)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: fileID
+      CHARACTER(*), INTENT(in):: title
+      INTEGER, INTENT(in):: iteration, dimensions, nNodes
+      REAL(8), INTENT(in):: time
+
+
+      WRITE(fileID, "(A)") '$NodeData'
+      WRITE(fileID, "(I10)") 1
+      WRITE(fileID, "(A1, A, A1)") '"' , title , '"'
+      WRITE(fileID, "(I10)") 1
+      WRITE(fileID, "(ES20.6E3)") time
+      WRITE(fileID, "(I10)") 3
+      WRITE(fileID, "(I10)") iteration
+      WRITE(fileID, "(I10)") dimensions
+      WRITE(fileID, "(I10)") nNodes
+
+    END SUBROUTINE writeGmsh2HeaderNodeData
+
+    !Footer
+    SUBROUTINE writeGmsh2FooterNodeData(fileID)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: fileID
+
+      WRITE(fileID, "(A)") '$EndNodeData'
+
+    END SUBROUTINE writeGmsh2FooterNodeData
+
+    !Element data subroutines
+    !Header
+    SUBROUTINE writeGmsh2HeaderElementData(fileID, title, iteration, time, dimensions, nVols)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: fileID
+      CHARACTER(*), INTENT(in):: title
+      INTEGER, INTENT(in):: iteration, dimensions, nVols
+      REAL(8), INTENT(in):: time
+
+
+      WRITE(fileID, "(A)") '$ElementData'
+      WRITE(fileID, "(I10)") 1
+      WRITE(fileID, "(A1, A, A1)") '"' , title , '"'
+      WRITE(fileID, "(I10)") 1
+      WRITE(fileID, "(ES20.6E3)") time
+      WRITE(fileID, "(I10)") 3
+      WRITE(fileID, "(I10)") iteration
+      WRITE(fileID, "(I10)") dimensions
+      WRITE(fileID, "(I10)") nVols
+
+    END SUBROUTINE writeGmsh2HeaderElementData
+
+    !Footer
+    SUBROUTINE writeGmsh2FooterElementData(fileID)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: fileID
+
+      WRITE(fileID, "(A)") '$EndElementData'
+
+    END SUBROUTINE writeGmsh2FooterElementData
+
     !Prints the scattered properties of particles into the nodes
     SUBROUTINE printOutputGmsh2(self, t)
       USE moduleMesh
@@ -21,65 +99,36 @@ MODULE moduleMeshOutputGmsh2
 
       DO i = 1, nSpecies
         WRITE(tstring, iterationFormat) t
-        fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
+        fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
-        WRITE(60, "(A)") '$MeshFormat'
+
-        WRITE(60, "(A)") '2.2 0 8'
+        CALL writeGmsh2HeaderMesh(60)
-        WRITE(60, "(A)") '$EndMeshFormat'
+
-        WRITE(60, "(A)") '$NodeData'
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes)
-        WRITE(60, "(A)") '1'
-        WRITE(60, "(A)") '"' // species(i)%obj%name // ' density (m^-3)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) self%numNodes
         DO n=1, self%numNodes
           CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
           WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
         END DO
-        WRITE(60, "(A)") '$EndNodeData'
+        CALL writeGmsh2FooterNodeData(60)
-        WRITE(60, "(A)") '$NodeData'
+
-        WRITE(60, "(A)") '1'
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes)
-        WRITE(60, "(A)") '"' // species(i)%obj%name // ' velocity (m s^-1)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 3
-        WRITE(60, *) self%numNodes
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
         END DO
-        WRITE(60, "(A)") '$EndNodeData'
+        CALL writeGmsh2FooterNodeData(60)
-        WRITE(60, "(A)") '$NodeData'
+
-        WRITE(60, "(A)") '1'
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes)
-        WRITE(60, "(A)") '"' // species(i)%obj%name // ' pressure (Pa)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) self%numNodes
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
         END DO
-        WRITE(60, "(A)") '$EndNodeData'
+        CALL writeGmsh2FooterNodeData(60)
-        WRITE(60, "(A)") '$NodeData'
+
-        WRITE(60, "(A)") '1'
+        CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes)
-        WRITE(60, "(A)") '"' // species(i)%obj%name // ' temperature (K)"'
-        WRITE(60, *) 1
-        WRITE(60, *) time
-        WRITE(60, *) 3
-        WRITE(60, *) t
-        WRITE(60, *) 1
-        WRITE(60, *) self%numNodes
         DO n=1, self%numNodes
           WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
         END DO
-        WRITE(60, "(A)") '$EndNodeData'
+        CALL writeGmsh2FooterNodeData(60)
         CLOSE (60)
 
       END DO
@@ -102,7 +151,9 @@ MODULE moduleMeshOutputGmsh2
       INTEGER:: n
       REAL(8):: time
       CHARACTER(:), ALLOCATABLE:: fileName
-      CHARACTER (LEN=iterationDigits):: tstring
+      CHARACTER (LEN=iterationDigits):: tString
+      CHARACTER(:), ALLOCATABLE:: title
+      CHARACTER (LEN=2):: cString
 
       SELECT TYPE(self)
       TYPE IS(meshParticles)
@@ -118,30 +169,26 @@ MODULE moduleMeshOutputGmsh2
 
       IF (collOutput) THEN
         time = DBLE(t)*tauMin*ti_ref
-        WRITE(tstring, iterationFormat) t
+        WRITE(tString, iterationFormat) t
 
-        fileName='OUTPUT_' // tstring// '_Collisions.msh'
+        fileName='OUTPUT_' // tString// '_Collisions.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (60, file = path // folder // '/' //  fileName)
-        WRITE(60, "(A)") '$MeshFormat'
+
-        WRITE(60, "(A)") '2.2 0 8'
+        CALL writeGmsh2HeaderMesh(60)
-        WRITE(60, "(A)") '$EndMeshFormat'
+
         DO k = 1, nCollPairs
           DO c = 1, interactionMatrix(k)%amount
-            WRITE(60, "(A)") '$ElementData'
+            WRITE(cString, "(I2)") c
-            WRITE(60, "(A)") '1'
+            title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
-            WRITE(60, "(5A,I2)") '"Pair ', interactionMatrix(k)%sp_i%name, '-', interactionMatrix(k)%sp_j%name, ' collision ', c
+            CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numVols)
-            WRITE(60, *) 1
-            WRITE(60, *) time
-            WRITE(60, *) 3
-            WRITE(60, *) t
-            WRITE(60, *) 1
-            WRITE(60, *) self%numVols
             DO n=1, self%numVols
               WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c)
             END DO
-            WRITE(60, "(A)") '$EndElementData'
+            CALL writeGmsh2FooterElementData(60)
+
           END DO
+
         END DO
 
         CLOSE(60)
@@ -175,50 +222,26 @@ MODULE moduleMeshOutputGmsh2
         fileName='OUTPUT_' // tstring// '_EMField.msh'
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (20, file = path // folder // '/' //  fileName)
-        WRITE(20, "(A)") '$MeshFormat'
+
-        WRITE(20, "(A)") '2.2 0 8'
+        CALL writeGmsh2HeaderMesh(20)
-        WRITE(20, "(A)") '$EndMeshFormat'
+
-        WRITE(20, "(A)") '$NodeData'
+        CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes)
-        WRITE(20, "(A)") '1'
-        WRITE(20, "(A)") '"Potential (V)"'
-        WRITE(20, *) 1
-        WRITE(20, *) time
-        WRITE(20, *) 3
-        WRITE(20, *) t
-        WRITE(20, *) 1
-        WRITE(20, *) self%numNodes
         DO n=1, self%numNodes
           WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
         END DO
-        WRITE(20, "(A)") '$EndNodeData'
+        CALL writeGmsh2FooterNodeData(20)
         
-        WRITE(20, "(A)") '$ElementData'
+        CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numVols)
-        WRITE(20, "(A)") '1'
-        WRITE(20, "(A)") '"Electric Field (V m^-1)"'
-        WRITE(20, *) 1
-        WRITE(20, *) time
-        WRITE(20, *) 3
-        WRITE(20, *) t
-        WRITE(20, *) 3
-        WRITE(20, *) self%numVols
         DO e=1, self%numVols
           WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
         END DO
-        WRITE(20, "(A)") '$EndElementData'
+        CALL writeGmsh2FooterElementData(20)
 
-        WRITE(20, "(A)") '$NodeData'
+        CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes)
-        WRITE(20, "(A)") '1'
-        WRITE(20, "(A)") '"Magnetic Field (T)"'
-        WRITE(20, *) 1
-        WRITE(20, *) time
-        WRITE(20, *) 3
-        WRITE(20, *) t
-        WRITE(20, *) 3
-        WRITE(20, *) self%numNodes
         DO n=1, self%numNodes
           WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref
         END DO
-        WRITE(20, "(A)") '$EndNodeData'
+        CALL writeGmsh2FooterNodeData(20)
 
         CLOSE(20)
 
@@ -226,4 +249,76 @@ MODULE moduleMeshOutputGmsh2
 
     END SUBROUTINE printEMGmsh2
 
+    !Prints the average properties of particles into the nodes
+    SUBROUTINE printAverageGmsh2(self)
+      USE moduleMesh
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleOutput
+      USE moduleAverage
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER:: n, i
+      TYPE(outputFormat), DIMENSION(1:self%numNodes):: outputMean, outputDeviation
+      CHARACTER(:), ALLOCATABLE:: fileName
+      INTEGER:: fileMean=10, fileDeviation=20
+
+      DO i = 1, nSpecies
+        fileName= 'Average_mean_' // species(i)%obj%name // '.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (fileMean, file = path // folder // '/' //  fileName)
+
+        fileName= 'Average_deviation_' // species(i)%obj%name // '.msh'
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (filedeviation, file = path // folder // '/' //  fileName)
+
+        CALL writeGmsh2HeaderMesh(fileMean)
+        CALL writeGmsh2HeaderMesh(fileDeviation)
+
+        CALL writeGmsh2HeaderNodeData(fileMean,      species(i)%obj%name // ' density, mean (m^-3)', 0, 0.D0, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' density, sd (m^-3)', 0, 0.D0, 1, self%numNodes)
+        DO n=1, self%numNodes
+          CALL calculateOutput(averageScheme(n)%mean%output(i),      outputMean(n),      self%nodes(n)%obj%v, species(i)%obj)
+          WRITE(fileMean,      "(I6,ES20.6E3)") n, outputMean(n)%density
+          CALL calculateOutput(averageScheme(n)%deviation%output(i), outputDeviation(n), self%nodes(n)%obj%v, species(i)%obj)
+          WRITE(fileDeviation, "(I6,ES20.6E3)") n, outputDeviation(n)%density
+        END DO
+        CALL writeGmsh2FooterNodeData(fileMean)
+        CALL writeGmsh2FooterNodeData(fileDeviation)
+
+        CALL writeGmsh2HeaderNodeData(fileMean,      species(i)%obj%name // ' velocity, mean (m s^-1)', 0, 0.D0, 3, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' velocity, sd (m s^-1)', 0, 0.D0, 3, self%numNodes)
+        DO n=1, self%numNodes
+          WRITE(fileMean,      "(I6,3(ES20.6E3))") n, outputMean(n)%velocity
+          WRITE(fileDeviation, "(I6,3(ES20.6E3))") n, outputDeviation(n)%velocity
+        END DO
+        CALL writeGmsh2FooterNodeData(fileMean)
+        CALL writeGmsh2FooterNodeData(fileDeviation)
+
+        CALL writeGmsh2HeaderNodeData(fileMean,      species(i)%obj%name // ' Pressure, mean (Pa)', 0, 0.D0, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' Pressure, sd (Pa)', 0, 0.D0, 1, self%numNodes)
+        DO n=1, self%numNodes
+          WRITE(fileMean,      "(I6,3(ES20.6E3))") n, outputMean(n)%pressure
+          WRITE(fileDeviation, "(I6,3(ES20.6E3))") n, outputDeviation(n)%pressure
+        END DO
+        CALL writeGmsh2FooterNodeData(fileMean)
+        CALL writeGmsh2FooterNodeData(fileDeviation)
+
+        CALL writeGmsh2HeaderNodeData(fileMean,      species(i)%obj%name // ' Temperature, mean (K)', 0, 0.D0, 1, self%numNodes)
+        CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' Temperature, sd (K)', 0, 0.D0, 1, self%numNodes)
+        DO n=1, self%numNodes
+          WRITE(fileMean,      "(I6,3(ES20.6E3))") n, outputMean(n)%temperature
+          WRITE(fileDeviation, "(I6,3(ES20.6E3))") n, outputDeviation(n)%temperature
+        END DO
+        CALL writeGmsh2FooterNodeData(fileMean)
+        CALL writeGmsh2FooterNodeData(fileDeviation)
+
+        CLOSE (fileMean)
+        CLOSE(fileDeviation)
+
+      END DO
+
+    END SUBROUTINE printAverageGmsh2
+
 END MODULE moduleMeshOutputGmsh2

src/modules/mesh/moduleMesh.f90

@@ -289,7 +289,7 @@ MODULE moduleMesh
     CONTAINS
       PROCEDURE, PASS:: doCollisions
 
-  END TYPE
+  END TYPE meshGeneric
 
   ABSTRACT INTERFACE
     !Reads the mesh from a file
@@ -338,9 +338,10 @@ MODULE moduleMesh
     REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
     !Permutation matrix for P L U factorization
     INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
-    PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
+    PROCEDURE(printOutput_interface),  POINTER, PASS:: printOutput  => NULL()
-    PROCEDURE(printEM_interface),     POINTER, PASS:: printEM     => NULL()
+    PROCEDURE(printEM_interface),      POINTER, PASS:: printEM      => NULL()
-    PROCEDURE(doCoulomb_interface),   POINTER, PASS:: doCoulomb   => NULL()
+    PROCEDURE(doCoulomb_interface),    POINTER, PASS:: doCoulomb    => NULL()
+    PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
     CONTAINS
       PROCEDURE, PASS:: constructGlobalK
 
@@ -373,6 +374,14 @@ MODULE moduleMesh
 
     END SUBROUTINE printEM_interface
 
+    !Prints average values
+    SUBROUTINE printAverage_interface(self)
+      IMPORT meshParticles
+
+      CLASS(meshParticles), INTENT(in):: self
+
+    END SUBROUTINE printAverage_interface
+
 
   END INTERFACE
 

src/modules/moduleAverage.f90

@@ -0,0 +1,72 @@
+MODULE moduleAverage
+  USE moduleOutput
+
+  TYPE:: averageData
+    TYPE(outputNode), ALLOCATABLE:: output(:)
+    TYPE(emNode):: emData
+
+  END TYPE averageData
+
+  !Generic type for average scheme
+  TYPE, PUBLIC:: averageGeneric
+    TYPE(averageData):: mean
+    TYPE(averageData):: deviation
+    CONTAINS
+      PROCEDURE, PASS:: firstAverage
+      PROCEDURE, PASS:: updateAverage
+
+  END TYPE averageGeneric
+
+  TYPE(averageGeneric), ALLOCATABLE:: averageScheme(:)
+
+  !Logical to determine if average scheme must be used
+  LOGICAL:: useAverage = .FALSE.
+  INTEGER:: tAverageStart = 1 !Starting iteartion for average scheme
+
+  CONTAINS
+    PURE SUBROUTINE firstAverage(self, output)
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(averageGeneric), INTENT(inout):: self
+      TYPE(outputNode), INTENT(in):: output(:)
+
+      !Copy current data as mean
+      self%mean%output(:) = output(:)
+
+      !Set deviation at 0
+      self%deviation%output(:) = 0.D0
+
+    END SUBROUTINE firstAverage
+
+
+    !Based on Welford's online algorithm
+    SUBROUTINE updateAverage(self, output, tAverage)
+      USE moduleOutput
+      USE moduleSpecies, ONLY: nSpecies
+      IMPLICIT NONE
+
+      CLASS(averageGeneric), INTENT(inout):: self
+      TYPE(outputNode), INTENT(in):: output(:)
+      INTEGER, INTENT(in):: tAverage
+      TYPE(averageData):: newMean, newDeviation
+
+      ALLOCATE(newMean%output(1:nSpecies))
+      ALLOCATE(newDeviation%output(1:nSpecies))
+
+      newMean%output(:) =  self%mean%output(:) + (output(:) - self%mean%output(:))/tAverage
+
+      newDeviation%output(:) = self%deviation%output(:) + ((output(:) - self%mean%output(:)) * &
+                                                           (output(:) - newMean%output(:))   - &
+                                                           self%deviation%output(:))/tAverage
+
+      self%mean%output(:) = newMean%output(:)
+
+      self%deviation%output(:) = newDeviation%output(:)
+
+      DEALLOCATE(newMean%output)
+      DEALLOCATE(newDeviation%output)
+        
+    END SUBROUTINE updateAverage
+
+END MODULE moduleAverage

src/modules/moduleInput.f90

@@ -75,6 +75,11 @@ MODULE moduleInput
       CALL readInject(config)
       CALL checkStatus(config, "readInject")
 
+      !Read average scheme
+      CALL verboseError('Reading average scheme...')
+      CALL readAverage(config)
+      CALL checkStatus(config, "readAverage")
+
       !Read parallel parameters
       CALL verboseError('Reading Parallel configuration...')
       CALL readParallel(config)
@@ -1180,6 +1185,40 @@ MODULE moduleInput
       
     END SUBROUTINE readInject
 
+    SUBROUTINE readAverage(config)
+      USE moduleAverage
+      USE moduleCaseParam, ONLY: tauMin
+      USE moduleMesh, ONLY: mesh
+      USE moduleSpecies, ONLY: nSpecies
+      IMPLICIT NONE
+
+      TYPE(json_file), INTENT(inout):: config
+      LOGICAL:: found
+      REAL(8):: tStart
+      INTEGER:: n
+
+      CALL config%info('average', found)
+      IF (found) THEN
+        useAverage = .TRUE.
+        CALL config%get('average.startTime', tStart, found)
+
+        IF (found) THEN
+          tAverageStart = INT(tStart / tauMin)
+
+        END IF
+
+        ALLOCATE(averageScheme(1:mesh%numNodes))
+
+        DO n = 1, mesh%numNodes
+          ALLOCATE(averageScheme(n)%mean%output(1:nSpecies))
+          ALLOCATE(averageScheme(n)%deviation%output(1:nSpecies))
+
+        END DO
+
+      END IF
+
+    END SUBROUTINE readAverage
+
     !Reads the velocity distribution functions for each inject
     SUBROUTINE readVelDistr(config, inj, object)
       USE moduleErrors

src/modules/moduleMath.f90

@@ -59,4 +59,15 @@ MODULE moduleMath
 
     END FUNCTION
 
+    FUNCTION tensorTrace(a) RESULT(t)
+      IMPLICIT NONE
+
+      REAL(8), DIMENSION(1:3,1:3):: a
+      REAL(8):: t
+
+      t = 0.D0
+      t = a(1,1)+a(2,2)+a(3,3)
+
+    END FUNCTION tensorTrace
+
 END MODULE moduleMath

src/modules/moduleOutput.f90

@@ -1,9 +1,22 @@
 !Contains information about output
 MODULE moduleOutput
   IMPLICIT NONE
+
   !Output for each node 
-  TYPE outputNode
+  TYPE, PUBLIC:: outputNode
     REAL(8):: den = 0.D0, mom(1:3) = 0.D0, tensorS(1:3,1:3) = 0.D0
+    CONTAINS
+      PROCEDURE, PASS(self), PRIVATE:: outputNode_equal_outputNode
+      PROCEDURE, PASS(self), PRIVATE:: outputNode_equal_real
+      PROCEDURE, PASS(self), PRIVATE:: outputNode_add_outputNode
+      PROCEDURE, PASS(self), PRIVATE:: outputNode_sub_outputNode
+      PROCEDURE, PASS(self), PRIVATE:: outputNode_mul_outputNode
+      PROCEDURE, PASS(self), PRIVATE:: outputNode_div_int
+      GENERIC, PUBLIC :: ASSIGNMENT(=) => outputNode_equal_outputNode, outputNode_equal_real
+      GENERIC, PUBLIC :: OPERATOR(+)   => outputNode_add_outputNode
+      GENERIC, PUBLIC :: OPERATOR(-)   => outputNode_sub_outputNode
+      GENERIC, PUBLIC :: OPERATOR(*)   => outputNode_mul_outputNode
+      GENERIC, PUBLIC :: OPERATOR(/)   => outputNode_div_int
 
   END TYPE
 
@@ -32,16 +45,81 @@ MODULE moduleOutput
   LOGICAL:: emOutput   = .FALSE.
 
   CONTAINS
-    FUNCTION tensorTrace(a) RESULT(t)
+    PURE SUBROUTINE outputNode_equal_outputNode(self, from)
       IMPLICIT NONE
 
-      REAL(8), DIMENSION(1:3,1:3):: a
+      CLASS(outputNode), INTENT(inout):: self
-      REAL(8):: t
+      CLASS(outputNode), INTENT(in):: from
 
-      t = 0.D0
+      self%den     = from%den
-      t = a(1,1)+a(2,2)+a(3,3)
+      self%mom     = from%mom
+      self%tensorS = from%tensorS
 
-    END FUNCTION tensorTrace
+    END SUBROUTINE outputNode_equal_outputNode
+
+    PURE ELEMENTAL SUBROUTINE outputNode_equal_real(self, from)
+      IMPLICIT NONE
+
+      CLASS(outputNode), INTENT(inout):: self
+      REAL(8), INTENT(in):: from
+
+      self%den     = from
+      self%mom     = from
+      self%tensorS = from
+
+    END SUBROUTINE outputNode_equal_real
+
+    PURE ELEMENTAL FUNCTION outputNode_add_outputNode(self, that) RESULT(total)
+      IMPLICIT NONE
+
+      CLASS(outputNode), INTENT(in):: self
+      CLASS(outputNode), INTENT(in):: that
+      TYPE(outputNode):: total
+
+      total%den     = self%den     + that%den
+      total%mom     = self%mom     + that%mom
+      total%tensorS = self%tensorS + that%tensorS
+
+    END FUNCTION outputNode_add_outputNode
+
+    PURE ELEMENTAL FUNCTION outputNode_sub_outputNode(self, that) RESULT(total)
+      IMPLICIT NONE
+
+      CLASS(outputNode), INTENT(in):: self
+      CLASS(outputNode), INTENT(in):: that
+      TYPE(outputNode):: total
+
+      total%den     = self%den     - that%den
+      total%mom     = self%mom     - that%mom
+      total%tensorS = self%tensorS - that%tensorS
+
+    END FUNCTION outputNode_sub_outputNode
+
+    PURE ELEMENTAL FUNCTION outputNode_mul_outputNode(self, that) RESULT(total)
+      IMPLICIT NONE
+
+      CLASS(outputNode), INTENT(in):: self
+      CLASS(outputNode), INTENT(in):: that
+      TYPE(outputNode):: total
+
+      total%den     = self%den     * that%den
+      total%mom     = self%mom     * that%mom
+      total%tensorS = self%tensorS * that%tensorS
+
+    END FUNCTION outputNode_mul_outputNode
+
+    PURE ELEMENTAL FUNCTION outputNode_div_int(self, that) RESULT(total)
+      IMPLICIT NONE
+
+      CLASS(outputNode), INTENT(in):: self
+      INTEGER, INTENT(in):: that
+      TYPE(outputNode):: total
+
+      total%den     = self%den     / REAL(that)
+      total%mom     = self%mom     / REAL(that)
+      total%tensorS = self%tensorS / REAL(that)
+
+    END FUNCTION outputNode_div_int
 
     SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn)
       USE moduleConstParam

src/modules/moduleSolver.f90

@@ -1,4 +1,5 @@
 MODULE moduleSolver
+  USE moduleAverage
 
   !Generic type for pusher of particles
   TYPE, PUBLIC:: pusherGeneric
@@ -15,6 +16,7 @@ MODULE moduleSolver
   !Generic type for solver
   TYPE, PUBLIC:: solverGeneric
     TYPE(pusherGeneric), ALLOCATABLE:: pusher(:)
+    TYPE(averageGeneric), ALLOCATABLE:: averageScheme
     PROCEDURE(solveEM_interface), POINTER, NOPASS:: solveEM => NULL()
     PROCEDURE(weightingScheme_interface), POINTER, NOPASS:: weightingScheme => NULL()
     CONTAINS
@@ -658,24 +660,18 @@ MODULE moduleSolver
       TYPE(particle), INTENT(inout):: part
       CLASS(meshVol), POINTER, INTENT(in):: volOld
       CLASS(meshVol), POINTER, INTENT(inout):: volNew
-      REAL(8):: fractionVolume, fractionWeight
+      REAL(8):: fractionVolume, pSplit
-      INTEGER:: nSplit
 
-      !If particle has change cell, call Weighting scheme
+      !If particle changes volume to smaller cell
-      IF (volOld%n /= volNew%n) THEN
+      IF (volOld%volume > volNew%volume) THEN
         fractionVolume = volOld%volume/volNew%volume
 
-        part%weight = part%weight * fractionVolume
+        !Calculate probability of splitting particle
+        pSplit = 1.D0 - DEXP(-fractionVolume)
 
-        fractionWeight = part%weight / part%species%weight
+        IF (random() < pSplit THEN
-
+          !Split particle in two
-        IF (fractionWeight >= 2.D0) THEN
+          CALL splitParticle(part, 2, volNew)
-          nSplit = FLOOR(fractionWeight)
-          CALL splitParticle(part, nSplit, volNew)
-
-        ELSEIF (part%weight < 1.D0) THEN
-          !Particle has lost statistical meaning and will be terminated
-          part%n_in = .FALSE.
 
         END IF
 
@@ -816,7 +812,44 @@ MODULE moduleSolver
 
       END IF
 
+      !Output average values
+      IF (useAverage .AND. t == tFinal) THEN
+        CALL mesh%printAverage()
+
+      END IF
+
     END SUBROUTINE doOutput
 
+    SUBROUTINE doAverage(t)
+      USE moduleAverage
+      USE moduleMesh
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: t
+      INTEGER:: tAverage, n
+
+
+      IF (useAverage) THEN
+        tAverage = t - tAverageStart
+
+        IF (tAverage == 1) THEN
+          !First iteration in which average scheme is used
+          DO n = 1, mesh%numNodes
+            CALL averageScheme(n)%firstAverage(mesh%nodes(n)%obj%output)
+
+          END DO
+        ELSEIF (tAverage > 1) THEN
+          DO n = 1, mesh%numNodes
+            CALL averageScheme(n)%updateAverage(mesh%nodes(n)%obj%output, tAverage)
+
+          END DO
+
+        END IF
+
+      END IF
+
+    END SUBROUTINE doAverage
+
+
 END MODULE moduleSolver