Modification to Probes to use all particles

Now, probes check all particles in the domain.
This is done only when probes are outputed to save CPU time.
However, still some issues and distribution functions are not properly
being calculated.

src/fpakc.f90

@@ -55,6 +55,9 @@ PROGRAM fpakc
 
     !Checks if a species needs to me moved in this iteration
     CALL solver%updatePushSpecies(t)
+
+    !Checks if probes need to be calculated this iteration
+    CALL resetProbes(t)
     tPush = omp_get_wtime()
     !$OMP END SINGLE
 
@@ -85,9 +88,6 @@ PROGRAM fpakc
     !$OMP SINGLE
     tCoul = omp_get_wTime() - tCoul
 
-    !Do probing
-    CALL doProbes(t)
-
     !Reset particles
     tReset = omp_get_wtime()
     !$OMP END SINGLE

src/modules/makefile

@@ -22,7 +22,7 @@ moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90
 moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleOutput.o: moduleMath.o moduleSpecies.o moduleRefParam.o moduleOutput.f90
+moduleOutput.o: moduleMath.o moduleSpecies.o moduleRefParam.o moduleProbe.o moduleOutput.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleRandom.o: moduleConstParam.o moduleRandom.f90
@@ -34,7 +34,7 @@ moduleRefParam.o: moduleConstParam.o moduleRefParam.f90
 moduleSpecies.o: moduleErrors.o moduleCaseParam.o moduleSpecies.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleSolver.o: mesh.o moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleOutput.o moduleSolver.f90
+moduleSolver.o: mesh.o moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleProbe.o moduleOutput.o moduleSolver.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleTable.o: moduleErrors.o moduleTable.f90

src/modules/moduleProbe.f90

@@ -1,11 +1,9 @@
 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
@@ -13,6 +11,9 @@ MODULE moduleProbe
     REAL(8), DIMENSION(1:3):: vrange
     REAL(8):: dvInv
     INTEGER:: every
+    LOGICAL:: update
+    INTEGER(KIND=OMP_LOCK_KIND):: lock
+    REAL(8):: maxR
     CONTAINS
       PROCEDURE, PASS:: init
       PROCEDURE, PASS:: findLowerIndex
@@ -41,7 +42,7 @@ MODULE moduleProbe
       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
+      INTEGER:: sp, i
       REAL(8):: dv(1:3)
 
       !Assign id
@@ -53,9 +54,6 @@ MODULE moduleProbe
 
       !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
@@ -95,6 +93,13 @@ MODULE moduleProbe
       !Number of iterations between output
       self%every = NINT(timeStep/ tauMin / ti_ref)
       
+      !Maximum radius
+      !TODO: Make this an input parameter
+      self%maxR = 1.D-1/L_ref
+
+      !Init the probe lock
+      CALL OMP_INIT_LOCK(self%lock)
+
     END SUBROUTINE init
 
     SUBROUTINE findLowerIndex(self, vp, i, j, k, inside)
@@ -114,32 +119,28 @@ MODULE moduleProbe
 
     END SUBROUTINE findLowerIndex
 
-    SUBROUTINE calculate(self)
+    SUBROUTINE calculate(self, part)
       USE moduleSpecies
       USE moduleList
       IMPLICIT NONE
 
       CLASS(probeDistFunc), INTENT(inout):: self
-      TYPE(particle), POINTER:: part
+      TYPE(particle), INTENT(in):: part
-      TYPE(lNode), POINTER:: node
+      REAL(8):: deltaR
       INTEGER:: i, j, k
       LOGICAL:: inside
+      REAL(8):: weight
       REAL(8):: fi, fi1
       REAL(8):: fj, fj1
       REAL(8):: fk, fk1
 
-      !Reset distribution function
+      !If particle is of desired type, include in the distribution function
-      self%f = 0.D0
+      IF (part%species%n == self%species%n) THEN
+        !Calculate distance between particle and probe
+        deltaR = NORM2(self%r - part%r)
 
-      !Loop over particles in cell
+        !Only include particle if it is inside the maximum radius
-      node => self%cell%listPart_in%head
+        IF (deltaR < self%maxR) THEN
-      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)
 
@@ -153,27 +154,30 @@ MODULE moduleProbe
             fk  = self%vk(k+1) - part%v(3)
             fk1 = part%v(3)    - self%vk(k)
 
+            ! weight = part%weight * DEXP(deltaR/self%maxR)
+            weight = part%weight
+
+            !Lock the probe
+            CALL OMP_SET_LOCK(self%lock)
+            
             !Assign particle weight to distribution function
-            self%f(i  , j  , k  ) = fi  * fj  * fk  * part%weight
+            self%f(i  , j  , k  ) = fi  * fj  * fk  * weight
-            self%f(i+1, j  , k  ) = fi1 * fj  * fk  * part%weight
+            self%f(i+1, j  , k  ) = fi1 * fj  * fk  * weight
-            self%f(i  , j+1, k  ) = fi  * fj1 * fk  * part%weight
+            self%f(i  , j+1, k  ) = fi  * fj1 * fk  * weight
-            self%f(i+1, j+1, k  ) = fi1 * fj1 * fk  * part%weight
+            self%f(i+1, j+1, k  ) = fi1 * fj1 * fk  * weight
-            self%f(i  , j  , k+1) = fi  * fj  * fk1 * part%weight
+            self%f(i  , j  , k+1) = fi  * fj  * fk1 * weight
-            self%f(i+1, j  , k+1) = fi1 * fj  * fk1 * part%weight
+            self%f(i+1, j  , k+1) = fi1 * fj  * fk1 * weight
-            self%f(i  , j+1, k+1) = fi  * fj1 * fk1 * part%weight
+            self%f(i  , j+1, k+1) = fi  * fj1 * fk1 * weight
-            self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * part%weight
+            self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight
+
+            !Unlock the probe
+            CALL OMP_UNSET_LOCK(self%lock)
 
           END IF
 
         END IF
 
-        !Move to next particle in the list
+      END IF
-        node => node%next
-
-      END DO
-
-      !Divide by the velocity cube volume
-      self%f = self%f * self%dvInv
           
     END SUBROUTINE calculate
 
@@ -182,13 +186,16 @@ MODULE moduleProbe
       USE moduleRefParam
       IMPLICIT NONE
 
-      CLASS(probeDistFunc), INTENT(in):: self
+      CLASS(probeDistFunc), INTENT(inout):: self
       INTEGER, INTENT(in):: t
       CHARACTER (LEN=iterationDigits):: tstring
       CHARACTER (LEN=3):: pstring
       CHARACTER(:), ALLOCATABLE:: filename
       INTEGER:: i, j, k
 
+      !Divide by the velocity cube volume
+      self%f = self%f * self%dvInv
+
       WRITE(tstring, iterationFormat) t
       WRITE(pstring, "(I3.3)") self%id
       fileName='OUTPUT_' // tstring// '_f_' // pstring // '.dat'
@@ -215,18 +222,21 @@ MODULE moduleProbe
 
       CLOSE(10)
 
+      !Reset distribution function
+      self%f = 0.D0
+
     END SUBROUTINE output
 
-    SUBROUTINE doProbes(t)
+    SUBROUTINE doProbes(part)
       IMPLICIT NONE
 
-      INTEGER, INTENT(in):: t
+      TYPE(particle), INTENT(in):: part
       INTEGER:: i
 
-      DO i = 1, SIZE(probe)
+      !Do it so they are only calculated when output
-        IF (MOD(t, probe(i)%every) == 0 .OR. t == tFinal) THEN
+      DO i = 1, nProbes
-          CALL probe(i)%calculate()
+        IF (probe(i)%update) THEN
-          CALL probe(i)%output(t)
+          CALL probe(i)%calculate(part)
 
         END IF
 
@@ -234,4 +244,34 @@ MODULE moduleProbe
 
     END SUBROUTINE doProbes
 
+    SUBROUTINE outputProbes(t)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: t
+      INTEGER:: i
+
+      DO i = 1, nProbes
+        IF (probe(i)%update) THEN
+          CALL probe(i)%output(t)
+
+        END IF
+
+      END DO
+
+    END SUBROUTINE outputProbes
+
+    SUBROUTINE resetProbes(t)
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: t
+      INTEGER:: i
+
+      DO i = 1, nProbes
+        probe(i)%f = 0.D0
+        probe(i)%update = MOD(t, probe(i)%every) == 0 .OR. t == tFinal
+
+      END DO
+
+    END SUBROUTINE resetProbes
+
 END MODULE moduleProbe

src/modules/moduleSolver.f90

@@ -497,17 +497,48 @@ MODULE moduleSolver
 
     END SUBROUTINE push0D
 
+    !Takes the particles from a list and put them into an array
+    !nStart indicates the last fill index in the array
+    SUBROUTINE resetList(partList, partArray, nStart)
+      USE moduleSpecies
+      USE moduleList
+      USE moduleProbe
+      IMPLICIT NONE
+
+      TYPE(listNode), INTENT(inout):: partList
+      TYPE(particle), INTENT(inout), DIMENSION(1:):: partArray
+      INTEGER, INTENT(in):: nStart
+      TYPE(lNode), POINTER:: partCurr, partNext
+      INTEGER:: n
+
+      partCurr => partList%head
+      DO n = 1, partList%amount
+       partNext => partCurr%next
+       partArray(nStart + n) = partCurr%part
+       CALL doProbes(partArray(nStart+n))
+       DEALLOCATE(partCurr)
+       partCurr => partNext
+
+      END DO
+
+      !Resest head, tail and number of nodes in the list
+      IF (ASSOCIATED(partList%head)) NULLIFY(partList%head)
+      IF (ASSOCIATED(partList%tail)) NULLIFY(partList%tail)
+      partList%amount = 0
+
+    END SUBROUTINE resetList
+
     SUBROUTINE doReset()
       USE moduleSpecies
       USE moduleMesh
       USE moduleList
+      USE moduleProbe
       IMPLICIT NONE
 
       INTEGER:: nn, n, e
       INTEGER, SAVE:: nPartNew
       INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces
       TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
-      TYPE(lNode), POINTER:: partCurr, partNext
 
       !$OMP SECTIONS
       !$OMP SECTION
@@ -524,10 +555,13 @@ MODULE moduleSolver
       END IF
       !$OMP SECTION
       nWScheme = partWScheme%amount
+
       !$OMP SECTION
       nCollisions = partCollisions%amount
+
       !$OMP SECTION
       nSurfaces = partSurfaces%amount
+
       !$OMP END SECTIONS
 
       !$OMP BARRIER
@@ -546,6 +580,7 @@ MODULE moduleSolver
         IF (partInj(n)%n_in) THEN
           nn = nn + 1
           partOld(nn) = partInj(n)
+          CALL doProbes(partOld(nn))
 
         END IF
 
@@ -558,54 +593,26 @@ MODULE moduleSolver
         IF (partTemp(n)%n_in) THEN
           nn = nn + 1
           partOld(nn) = partTemp(n)
+          CALL doProbes(partOld(nn))
 
         END IF
 
       END DO
+
       !$OMP SECTION
       !Reset particles from weighting scheme
       nn = nInjIn + nOldIn
-      partCurr => partWScheme%head
+      CALL resetList(partWScheme, partOld, nn)
-      DO n = 1, nWScheme
-        partNext => partCurr%next
-        partOld(nn+n) = partCurr%part
-        DEALLOCATE(partCurr)
-        partCurr => partNext
-
-      END DO
-      IF (ASSOCIATED(partWScheme%head)) NULLIFY(partWScheme%head)
-      IF (ASSOCIATED(partWScheme%tail)) NULLIFY(partWScheme%tail)
-      partWScheme%amount = 0
 
       !$OMP SECTION
       !Reset particles from collisional process
       nn = nInjIn + nOldIn + nWScheme
-      partCurr => partCollisions%head
+      CALL resetList(partCollisions, partOld, nn)
-      DO n = 1, nCollisions
-        partNext => partCurr%next
-        partOld(nn+n) = partCurr%part
-        DEALLOCATE(partCurr)
-        partCurr => partNext
-
-      END DO
-      IF (ASSOCIATED(partCollisions%head)) NULLIFY(partCollisions%head)
-      IF (ASSOCIATED(partCollisions%tail)) NULLIFY(partCollisions%tail)
-      partCollisions%amount = 0
 
       !$OMP SECTION
       !Reset particles from surface process
       nn = nInjIn + nOldIn + nWScheme + nCollisions
-      partCurr => partSurfaces%head
+      CALL resetList(partSurfaces, partOld, nn)
-      DO n = 1, nSurfaces
-        partNext => partCurr%next
-        partOld(nn+n) = partCurr%part
-        DEALLOCATE(partCurr)
-        partCurr => partNext
-
-      END DO
-      IF (ASSOCIATED(partSurfaces%head)) NULLIFY(partSurfaces%head)
-      IF (ASSOCIATED(partSurfaces%tail)) NULLIFY(partSurfaces%tail)
-      partSurfaces%amount = 0
 
       !$OMP SECTION
       !Reset output in nodes
@@ -827,6 +834,7 @@ MODULE moduleSolver
       USE moduleOutput
       USE moduleSpecies
       USE moduleCompTime
+      USE moduleProbe
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: t
@@ -838,6 +846,7 @@ MODULE moduleSolver
         !Resets output counter
         counterOutput=0
 
+        CALL outputProbes(t)
         CALL mesh%printOutput(t)
         IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
         CALL mesh%printEM(t)