From e683c66ff899da905e0e6d8b1664ecb0ea566d26 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 18 Oct 2021 16:00:28 +0200 Subject: [PATCH] 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 | 6 +- src/modules/makefile | 4 +- src/modules/moduleProbe.f90 | 128 +++++++++++++++++++++++------------ src/modules/moduleSolver.f90 | 77 +++++++++++---------- 4 files changed, 132 insertions(+), 83 deletions(-) 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)