diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 8d33917..03d0bbb 100644 --- a/src/fpakc.f90 +++ b/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 diff --git a/src/modules/makefile b/src/modules/makefile index 7f54dae..4eda2be 100644 --- a/src/modules/makefile +++ b/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 diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90 index 6093ac6..2d2295a 100644 --- a/src/modules/moduleProbe.f90 +++ b/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(lNode), POINTER:: node + TYPE(particle), INTENT(in):: part + 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 - self%f = 0.D0 + !If particle is of desired type, include in the distribution function + 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 - 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 + !Only include particle if it is inside the maximum radius + IF (deltaR < self%maxR) 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+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 + self%f(i , j , k ) = fi * fj * fk * weight + self%f(i+1, j , k ) = fi1 * fj * fk * weight + self%f(i , j+1, k ) = fi * fj1 * fk * weight + self%f(i+1, j+1, k ) = fi1 * fj1 * fk * weight + self%f(i , j , k+1) = fi * fj * fk1 * weight + self%f(i+1, j , k+1) = fi1 * fj * fk1 * weight + self%f(i , j+1, k+1) = fi * fj1 * fk1 * 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 - node => node%next - - END DO - - !Divide by the velocity cube volume - self%f = self%f * self%dvInv - + END IF + 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) - IF (MOD(t, probe(i)%every) == 0 .OR. t == tFinal) THEN - CALL probe(i)%calculate() - CALL probe(i)%output(t) + !Do it so they are only calculated when output + DO i = 1, nProbes + IF (probe(i)%update) THEN + 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 diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 0f78a5e..47d33ca 100644 --- a/src/modules/moduleSolver.f90 +++ b/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 - 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 + CALL resetList(partWScheme, partOld, nn) !$OMP SECTION !Reset particles from collisional process nn = nInjIn + nOldIn + nWScheme - partCurr => partCollisions%head - 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 + CALL resetList(partCollisions, partOld, nn) !$OMP SECTION !Reset particles from surface process nn = nInjIn + nOldIn + nWScheme + nCollisions - partCurr => partSurfaces%head - 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 + CALL resetList(partSurfaces, partOld, nn) !$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)