Merge branch 'feature/probe' into 'development'

Modification to Probes to use all particles

See merge request JorgeGonz/fpakc!25

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
@@ -94,6 +92,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
 
@@ -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,28 +154,31 @@ 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
 
     SUBROUTINE output(self, t)
@@ -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)