Merge branch 'feature/probe' into 'development'

First implementation of probing method

See merge request JorgeGonz/fpakc!23

doc/user-manual/fpakc_UserManual.tex

@@ -215,11 +215,12 @@
 
           \item Recombination.
                 When an electron and an ion interact, there is a possibility for them to be recombined into a neutral particle.
-                The photon emitted by this process is not modeled yet.
+                The photon emitted by this process is not modelled yet.
         \end{itemize}
       
-      \subsection{\acrlong{cs}}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-        Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline.
+%      \subsection{\acrlong{cs}}
+%        Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline.
         
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Reset of particle array}
@@ -228,6 +229,20 @@
       In this section, particles are assigned to the list of particles inside each individual cell.
       Unfortunately, this is done right now without parallelisation and is very CPU consuming.
       
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    \section{Probing}\label{sec:probing}
+      As default, \gls{fpakc} outputs information of macroscopic quantities (density, velocity, temperature\ldots) in the finite element mesh.
+      However, a lot of information can be extracted from the particle distribution function.
+      Thus, a probing method is provided to extract the distribution function in a specific position.
+      
+      The particles inside a cell in which the input position is located are distributed into a 3D velocity grid.
+      The user can decide the grid width and the number of points in each direction.
+      The distribution function will be calculated and wrote with a time step decided by the user.
+      
+      If a particle velocity resides outside of the velocity grid (in any direction), it wont be added to the tally of the distribution function.
+      Due to the limitation of only taking into account particles in the cell, and not neighbour particles, two probes for the same species at different positions but in the same cell will output the same results.
+      A more advance method taking into account distance between the particles and the probe position as well as particles in neighbour cells could be implemented to improve the statistics of the distribution function.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Scattering}
       The properties of each particle are deposited in the nodes from the containing cell.
@@ -359,6 +374,30 @@ make
                                   Trigger between writings is the same as in \textbf{triggerOutput}.
           \item \textbf{EMField}: Logical.
                                   Determines if the electromagnetic field is printed.
+          \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.
+                                 The object is structured as follows:
+                                 \begin{itemize}
+                                   \item \textbf{species}: Character.
+                                                           Species name as defined in \textbf{species} array.
+                                   \item \textbf{position}: Real.
+                                                            Array of dimension $3$.
+                                                            Units in $\unit{m}$.
+                                                            Indicates the position of the probing.
+                                   \item \textbf{timeStep}: Real
+                                                            Units in $\unit{s}$.
+                                                            Optional.
+                                                            Time step for output of the distribution function.
+                                                            If none is provided, the minimum time step of the case is used.
+                                   \item \textbf{velocity\_1, velocity\_2, velocity\_3}: Real.
+                                                                                         Array of dimension $2$.
+                                                                                         Velocity range (minimum-maximum) in which the distribution function will be interpolated.
+                                                                                         The subscripts 1, 2, 3 indicate the three directions of the case.
+                                   \item \textbf{points}: Integer.
+                                                          Array of dimension $3$.
+                                                          Number of points in each direction.
+                                 \end{itemize}
         \end{itemize}
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -439,9 +478,11 @@ make
         \begin{itemize}
           \item \textbf{name}: Character.
                                Name of the boundary.
-          \item \textbf{type}: Character.
+          \item \textbf{physicalSurface}: Integer.
-                               Type of boundary.
+                                          Identification of the surface in the mesh file.
-                               Accepted values are:
+          \item \textbf{bType}: Array of objects of dimension 'number of species'.
+                                Per each species defined in the case, a boundary \textbf{type} needs to be provided.
+                                Accepted values for \textbf{type} are:
                                 \begin{itemize}
                                   \item \textbf{reflection}: Elastic reflection of particles.
                                   \item \textbf{absorption}: Particle is eliminated from the domain.
@@ -503,8 +544,6 @@ make
                                   \item \textbf{axis}: Identifies the symmetry axis for 2D cylindrical simulations.
                                                        If for some reason a particle interact with this axis, it is reflected.
                                 \end{itemize}
-          \item \textbf{physicalSurface}: Integer.
-                                          Identification of the edge in the mesh file.
         \end{itemize}
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -631,8 +670,7 @@ make
                                    \begin{itemize}
                                       \item \textbf{Electrostatic}: Solves the Poison equation to obtain the self-consistent electrostatic potential.
                                    \end{itemize}
-          \item \textbf{initial}: Object.
+          \item \textbf{initial}: Array of objects.
-                                  Array.
                                   Determines initial values for the species.
                                   Required values are:
                                   \begin{itemize}
@@ -661,8 +699,7 @@ make
                                    Units of $\unit{s}$.
                                    Time step to calculate MCC.
                                    If no time is provided, the minimum time step is used.
-          \item \textbf{collisions}: Object.
+          \item \textbf{collisions}: Array of objects.
-                                     Array.
                                      Contains the different short range interactions (\acrshort{mcc}).
                                      Multiple collision types can be defined for each pair of species.
                                      Each object in the array is defined by:
@@ -670,8 +707,7 @@ make
                                        \item \textbf{species\_i}, \textbf{species\_j}: Character.
                                                                                        Define the two species involved in the collision processes.
                                                                                        Order is indiferent.
-                                       \item \textbf{cTypes}: Object.
+                                       \item \textbf{cTypes}: Array of objects.
-                                                              Array.
                                                               Defines all the collisions between \textbf{species\_i} and \textbf{species\_j}.
                                                               Required values are:
                                                               \begin{itemize}
@@ -705,29 +741,29 @@ make
                                  
         \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  \chapter{Example runs\label{ch:exampleRuns}}
+  \chapter{Example runs}\label{ch:exampleRuns}
     This chapter presents a description of the different example files distributed with \acrshort{fpakc}.
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{1D Emissive Cathode (1D\_Cathode)}
-      Emission from a 1D cathond in both, cartesian and radial coordinates.
+      Emission from a 1D cathode in both, cartesian and radial coordinates.
       Both cases insert the same amount of electrons from the minimum coordinate and have the same boundary conditions for particles and the electrostatic field.
-      This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifiying some parameters in the input file.
+      This case is useful to ilustrate hoy \acrshort{fpakc} can deal with different geometries by just modifying some parameters in the input file.
       The same mesh file (\lstinline|mesh.msh|) is used for both cases.
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{0D \ce{Ar}-\ce{Ar+} Elastic Collision (0D\_Argon)}
       Test case to check the 0D geometry that includes the elastic collision between \ce{Ar} and \ce{Ar+}.
-      Initial states are readed from the Argon\_Initial.dat and Argon+\_Initial.dat
+      Initial states are read from the Argon\_Initial.dat and Argon+\_Initial.dat
     
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{ALPHIE Grid system (ALPHIE\_Grid)}
-      Two-dimensional axialsymmetry case to study the counterflow of electrons and Argon ions going trhough the ALPHIE grid system.
+      Two-dimensional axial-symmetry case to study the counterflow of electrons and Argon ions going through the ALPHIE grid system.
       A \lstinline|mesh.geo| file is provided to easily modify the parameters of the grid system and generate a new mesh with \Gls{gmsh}.
       
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     \section{Flow around cylinder (cylFlow)}
-      Simple case of neutral Argon flow around a cylinder in a 2D axialsymmetry geometry.
+      Simple case of neutral Argon flow around a cylinder in a 2D axial-symmetry geometry.
       Elastic collisions between argon particles are included.
       Two cases are presented here: one in which the same mesh (meshSingle.msh) for scattering and collisions is used (input.json) and a second one (inputDualMesh.json) in which a mesh is used for scattering (mesh.msh) and a second one is used only for collisions (meshColl.msh).
     

src/fpakc.f90

@@ -1,12 +1,13 @@
 !  FPAKC main program
 PROGRAM fpakc
+  USE moduleCompTime
+  USE moduleCaseParam
   USE moduleInput
-  USE moduleErrors
   USE moduleInject
   USE moduleSolver
   USE moduleMesh
-  USE moduleCompTime
+  USE moduleProbe
-  USE moduleCaseParam
+  USE moduleErrors
   USE OMP_LIB
   IMPLICIT NONE
 
@@ -84,6 +85,9 @@ PROGRAM fpakc
     !$OMP SINGLE
     tCoul = omp_get_wTime() - tCoul
 
+    !Do probing
+    CALL doProbes(t)
+
     !Reset particles
     tReset = omp_get_wtime()
     !$OMP END SINGLE

src/makefile

@@ -4,6 +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)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
 			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \

src/modules/makefile

@@ -2,7 +2,8 @@
 OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
 		   moduleOutput.o moduleInput.o moduleSolver.o \
 			 moduleCollisions.o moduleTable.o moduleParallel.o \
-			 moduleEM.o moduleRandom.o moduleMath.o
+			 moduleEM.o moduleRandom.o moduleMath.o \
+			 moduleProbe.o
 
 all: $(OBJS) mesh.o
 
@@ -12,7 +13,7 @@ mesh.o: moduleCollisions.o moduleBoundary.o moduleMath.o
 moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleInput.f90
+moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleProbe.o moduleInput.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90

src/modules/moduleInput.f90

@@ -60,6 +60,11 @@ MODULE moduleInput
       CALL readGeometry(config)
       CALL checkStatus(config, "readGeometry")
 
+      !Read probes
+      CALL verboseError('Reading Probes...')
+      CALL readProbes(config)
+      CALL checkStatus(config, 'readProbes')
+
       !Read initial state for species
       CALL verboseError('Reading Initial state...')
       CALL readInitial(config)
@@ -903,6 +908,47 @@ MODULE moduleInput
 
     END SUBROUTINE readGeometry
 
+    SUBROUTINE readProbes(config)
+      USE moduleProbe
+      USE moduleErrors
+      USE json_module
+      IMPLICIT NONE
+
+      TYPE(json_file), INTENT(inout):: config
+      CHARACTER(:), ALLOCATABLE:: object
+      LOGICAL:: found
+      CHARACTER(2):: istring
+      INTEGER:: i
+      CHARACTER(:), ALLOCATABLE:: speciesName
+      REAL(8), ALLOCATABLE, DIMENSION(:):: r
+      REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
+      INTEGER, ALLOCATABLE, DIMENSION(:):: points
+      REAL(8):: timeStep
+
+      CALL config%info('output.probes', found, n_children = nProbes)
+
+      IF (found) ALLOCATE(probe(1:nProbes))
+
+      DO i = 1, nProbes
+        WRITE(iString, '(I2)') i
+        object = 'output.probes(' // trim(istring) // ')'
+
+        CALL config%get(object // '.species', speciesName, found)
+        CALL config%get(object // '.position', r, found)
+        CALL config%get(object // '.velocity_1', v1, found)
+        CALL config%get(object // '.velocity_2', v2, found)
+        CALL config%get(object // '.velocity_3', v3, found)
+        CALL config%get(object // '.points', points, found)
+        CALL config%get(object // '.timeStep', timeStep, found)
+
+        IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')
+
+        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep)
+
+      END DO
+
+    END SUBROUTINE readProbes
+
     SUBROUTINE readEMBoundary(config)
       USE moduleMesh
       USE moduleOutput
@@ -926,7 +972,7 @@ MODULE moduleInput
       IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
       
       DO i = 1, nBoundaryEM
-        WRITE(istring, '(i2)') i
+        WRITE(istring, '(I2)') i
         object = 'boundaryEM(' // trim(istring) // ')'
 
         CALL config%get(object // '.type', boundEM(i)%typeEM, found)

src/modules/moduleProbe.f90

@@ -0,0 +1,239 @@
+MODULE moduleProbe
+  USE moduleMesh
+  USE moduleSpecies
+
+  TYPE:: probeDistFunc
+    INTEGER:: id
+    REAL(8), DIMENSION(1:3):: r
+    CLASS(meshVol), POINTER:: cell
+    CLASS(speciesGeneric), POINTER:: species
+    REAL(8), ALLOCATABLE, DIMENSION(:):: vi, vj, vk
+    REAL(8), ALLOCATABLE, DIMENSION(:,:,:):: f
+    INTEGER, DIMENSION(1:3):: nv
+    REAL(8), DIMENSION(1:3):: vrange
+    REAL(8):: dvInv
+    INTEGER:: every
+    CONTAINS
+      PROCEDURE, PASS:: init
+      PROCEDURE, PASS:: findLowerIndex
+      PROCEDURE, PASS:: calculate
+      PROCEDURE, PASS:: output
+      
+  END TYPE probeDistFunc
+
+  INTEGER:: nProbes = 0
+  TYPE(probeDistFunc), ALLOCATABLE:: probe(:)
+
+  CONTAINS
+    !Functions for probeDistFunc type
+    SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep)
+      USE moduleCaseParam
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleMesh
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(probeDistFunc), INTENT(out):: self
+      INTEGER, INTENT(in):: id
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: speciesName
+      REAL(8), INTENT(in):: r(1:3)
+      REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2)
+      INTEGER, INTENT(in):: points(1:3)
+      REAL(8), INTENT(in):: timeStep
+      INTEGER:: sp, e, i
+      REAL(8):: dv(1:3)
+
+      !Assign id
+      self%id = id
+
+      !Get species
+      sp = speciesName2Index(speciesName)
+      self%species => species(sp)%obj
+
+      !Find cell
+      self%r = r/L_ref
+      e = findCellBrute(mesh, self%r)
+      IF (e == 0) CALL criticalError("No cell found for position in probe", 'init')
+      self%cell => mesh%vols(e)%obj
+
+      !Allocates velocity grid
+      self%nv = points
+      ALLOCATE(self%vi(1:self%nv(1)))
+      ALLOCATE(self%vj(1:self%nv(2)))
+      ALLOCATE(self%vk(1:self%nv(3)))
+
+      !Creates grid
+      dv(1) = (v1(2) - v1(1))/REAL(self%nv(1) - 1)/v_ref
+      dv(2) = (v2(2) - v2(1))/REAL(self%nv(2) - 1)/v_ref
+      dv(3) = (v3(2) - v3(1))/REAL(self%nv(3) - 1)/v_ref
+
+      self%dvInv = 1.D0/(dv(1)*dv(2)*dv(3))
+
+      DO i = 1, self%nv(1)
+        self%vi(i) = dv(1)*REAL(i - 1) + v1(1)/v_ref
+
+      END DO
+
+      DO i = 1, self%nv(2)
+        self%vj(i) = dv(2)*REAL(i - 1) + v2(1)/v_ref
+
+      END DO
+
+      DO i = 1, self%nv(3)
+        self%vk(i) = dv(3)*REAL(i - 1) + v3(1)/v_ref
+
+      END DO
+
+      self%vrange(1) = self%vi(self%nv(1)) - self%vi(1)
+      self%vrange(2) = self%vj(self%nv(2)) - self%vj(1)
+      self%vrange(3) = self%vk(self%nv(3)) - self%vk(1)
+
+      !Allocates distribution function
+      ALLOCATE(self%f(1:self%nv(1), &
+                      1:self%nv(2), &
+                      1:self%nv(3)))
+
+      !Number of iterations between output
+      self%every = NINT(timeStep/ tauMin / ti_ref)
+
+    END SUBROUTINE init
+
+    SUBROUTINE findLowerIndex(self, vp, i, j, k, inside)
+      IMPLICIT NONE
+      
+      CLASS(probeDistFunc), INTENT(in):: self
+      REAL(8), INTENT(in):: vp(1:3)
+      INTEGER, INTENT(out):: i, j, k
+      LOGICAL, INTENT(out):: inside
+
+      i = FLOOR((vp(1) - self%vi(1))/self%vrange(1)*(REAL(self%nv(1) - 1)) + 1.D0)
+      IF (i >= self%nv(1) .OR. i <  1) inside = .FALSE.
+      j = FLOOR((vp(2) - self%vj(1))/self%vrange(2)*(REAL(self%nv(2) - 1)) + 1.D0)
+      IF (j >= self%nv(2) .OR. j <  1) inside = .FALSE.
+      k = FLOOR((vp(3) - self%vk(1))/self%vrange(3)*(REAL(self%nv(3) - 1)) + 1.D0)
+      IF (k >= self%nv(3) .OR. k <  1) inside = .FALSE.
+
+    END SUBROUTINE findLowerIndex
+
+    SUBROUTINE calculate(self)
+      USE moduleSpecies
+      USE moduleList
+      IMPLICIT NONE
+
+      CLASS(probeDistFunc), INTENT(inout):: self
+      TYPE(particle), POINTER:: part
+      TYPE(lNode), POINTER:: node
+      INTEGER:: i, j, k
+      LOGICAL:: inside
+      REAL(8):: fi, fi1
+      REAL(8):: fj, fj1
+      REAL(8):: fk, fk1
+
+      !Reset distribution function
+      self%f = 0.D0
+
+      !Loop over particles in cell
+      node => self%cell%listPart_in%head
+      DO WHILE(ASSOCIATED(node))
+
+        !Selects particle on list
+        part => node%part
+
+        !If particle is of desired type, include in the distribution function
+        IF (part%species%n == self%species%n) THEN
+          !find lower index for all dimensions
+          CALL self%findLowerIndex(part%v, i, j, k, inside)
+
+          !If particle is inside the velocity grid, add it to the distribution function
+          IF (inside) THEN
+            !Calculate weights
+            fi  = self%vi(i+1) - part%v(1)
+            fi1 = part%v(1)    - self%vi(i)
+            fj  = self%vj(j+1) - part%v(2)
+            fj1 = part%v(2)    - self%vj(j)
+            fk  = self%vk(k+1) - part%v(3)
+            fk1 = part%v(3)    - self%vk(k)
+
+            !Assign particle weight to distribution function
+            self%f(i  , j  , k  ) = fi  * fj  * fk  * part%weight
+            self%f(i+1, j  , k  ) = fi1 * fj  * fk  * part%weight
+            self%f(i  , j+1, k  ) = fi  * fj1 * fk  * part%weight
+            self%f(i+1, j+1, k  ) = fi1 * fj1 * fk  * part%weight
+            self%f(i  , j  , k+1) = fi  * fj  * fk1 * part%weight
+            self%f(i+1, j  , k+1) = fi1 * fj  * fk1 * part%weight
+            self%f(i  , j+1, k+1) = fi  * fj1 * fk1 * part%weight
+            self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * part%weight
+
+          END IF
+          
+        END IF
+
+        !Move to next particle in the list
+        node => node%next
+
+      END DO
+
+      !Divide by the velocity cube volume
+      self%f = self%f * self%dvInv
+
+    END SUBROUTINE calculate
+
+    SUBROUTINE output(self, t)
+      USE moduleOutput
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(probeDistFunc), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      CHARACTER (LEN=iterationDigits):: tstring
+      CHARACTER (LEN=3):: pstring
+      CHARACTER(:), ALLOCATABLE:: filename
+      INTEGER:: i, j, k
+
+      WRITE(tstring, iterationFormat) t
+      WRITE(pstring, "(I3.3)") self%id
+      fileName='OUTPUT_' // tstring// '_f_' // pstring // '.dat'
+      WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+      OPEN (10, file = path // folder // '/' //  fileName)
+      WRITE(10, "(A1, 1X, A)")               "# ", self%species%name
+      WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m"
+      WRITE(10, "(A6, ES15.6E3, A2)")    "# t = ", REAL(t)*tauMin*ti_ref, " s"
+      WRITE(10, "(A1, A19, 3(A20))")     "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f"
+      DO i = 1, self%nv(1)
+        DO j = 1, self%nv(2)
+          DO k = 1, self%nv(3)
+            WRITE(10, "(4(ES20.6E3))") self%vi(i)*v_ref, &
+                                       self%vj(j)*v_ref, &
+                                       self%vk(k)*v_ref, &
+                                       self%f(i, j, k)
+
+          END DO
+          WRITE(10, *)
+
+        END DO
+
+      END DO
+
+      CLOSE(10)
+
+    END SUBROUTINE output
+
+    SUBROUTINE doProbes(t)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: t
+      INTEGER:: i
+
+      DO i = 1, SIZE(probe)
+        IF (MOD(t, probe(i)%every) == 0 .OR. t == tFinal) THEN
+          CALL probe(i)%calculate()
+          CALL probe(i)%output(t)
+
+        END IF
+
+      END DO
+
+    END SUBROUTINE doProbes
+
+END MODULE moduleProbe