From ffb03e634bd476add312ad3e23e1742cb7df39ab Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 13 Oct 2020 18:16:18 +0200 Subject: [PATCH] Minor improvements in performance and code clarity. Still no solution for the reset subroutine. It is really time consumming. --- makefile | 2 +- runs/cylFlow/compCPUTime.gp | 116 ++++++++++++++++++++++++------- runs/cylFlow/plotCPUTime.gp | 4 +- src/DSMC_Neutrals.f95 | 10 +-- src/modules/makefile | 4 +- src/modules/moduleCollisions.f95 | 14 ++-- src/modules/moduleInject.f95 | 2 + src/modules/moduleList.f95 | 32 +++++---- src/modules/moduleMesh.f95 | 17 +++-- src/modules/moduleMeshCyl.f95 | 79 ++++++++++++++------- src/modules/moduleSolver.f95 | 69 ++++++++++-------- src/modules/moduleSpecies.f95 | 1 - 12 files changed, 233 insertions(+), 117 deletions(-) diff --git a/makefile b/makefile index 749dad2..b10a4ac 100644 --- a/makefile +++ b/makefile @@ -13,7 +13,7 @@ INCDIR := INCDIR += $(JSONINC) # compile flags -FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall +FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -fno-stack-arrays -g #Output file OUTPUT = DSMC_Neutrals diff --git a/runs/cylFlow/compCPUTime.gp b/runs/cylFlow/compCPUTime.gp index 7b1155d..66a8bed 100644 --- a/runs/cylFlow/compCPUTime.gp +++ b/runs/cylFlow/compCPUTime.gp @@ -1,5 +1,5 @@ -set terminal qt enhanced 1 size 1600, 1000 font "Times ,10" +set terminal qt enhanced 1 persist size 1600, 1000 font "Times ,10" set style line 1 pt 4 lc rgb "#B50427" #Squares red set style line 2 pt 6 lc rgb "#3B4CC1" #Circles blue @@ -8,47 +8,109 @@ set style line 4 pt 2 lc rgb "#FE7F0E" #Exes orange set style line 5 pt 8 lc rgb "#D6696B" #Triangles light red set style line 10 lt 1 lw 2.0 lc rgb "black" #Black solid line -folder1 = "2020-09-27_10.46.31/" -folder2 = "2020-09-28_10.49.22/" -folder3 = "2020-09-30_20.40.04/" -folder4 = "2020-10-01_08.53.50/" +folder1 = "base_case/" +name1 = "Base Case" +folder2 = "2020-10-11_17.58.30/" +name2 = "New Version" +folder3 = "2020-10-13_17.53.32/" +name3 = "Sections" +folder4 = "2020-10-11_15.20.09/" +name4 = "New list" set key box opaque set pointsize 1.5 set ylabel "Time per particle (ms)" -set xrange [0:8000] +set xrange [0:500] set xlabel "Iteration" set multiplot layout 2,3 set title "Total" -plot folder1."cpuTime.dat" u 1:(1e3*$3) t "Step 1" ls 1, \ - folder2."cpuTime.dat" u 1:(1e3*$3) t "Step 2" ls 2, \ - folder3."cpuTime.dat" u 1:(1e3*$3) t "Step 3" ls 3, \ - folder4."cpuTime.dat" u 1:(1e3*$3) t "Step 4" ls 4 +plot folder1."cpuTime.dat" u 1:(1e3*$3) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e3*$3) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e3*$3) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e3*$3) t name4 ls 4 set title "Push" -plot folder1."cpuTime.dat" u 1:(1e3*$4) t "Step 1" ls 1, \ - folder2."cpuTime.dat" u 1:(1e3*$4) t "Step 2" ls 2, \ - folder3."cpuTime.dat" u 1:(1e3*$4) t "Step 3" ls 3, \ - folder4."cpuTime.dat" u 1:(1e3*$4) t "Step 4" ls 4 +plot folder1."cpuTime.dat" u 1:(1e3*$4) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e3*$4) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e3*$4) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e3*$4) t name4 ls 4 set title "Reset" -plot folder1."cpuTime.dat" u 1:(1e3*$5) t "Step 1" ls 1, \ - folder2."cpuTime.dat" u 1:(1e3*$5) t "Step 2" ls 2, \ - folder3."cpuTime.dat" u 1:(1e3*$5) t "Step 3" ls 3, \ - folder4."cpuTime.dat" u 1:(1e3*$5) t "Step 4" ls 4 +plot folder1."cpuTime.dat" u 1:(1e3*$5) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e3*$5) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e3*$5) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e3*$5) t name4 ls 4 set title "Collisions" -plot folder1."cpuTime.dat" u 1:(1e3*$6) t "Step 1" ls 1, \ - folder2."cpuTime.dat" u 1:(1e3*$6) t "Step 2" ls 2, \ - folder3."cpuTime.dat" u 1:(1e3*$6) t "Step 3" ls 3, \ - folder4."cpuTime.dat" u 1:(1e3*$6) t "Step 4" ls 4 +plot folder1."cpuTime.dat" u 1:(1e3*$6) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e3*$6) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e3*$6) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e3*$6) t name4 ls 4 set title "Weighting" -plot folder1."cpuTime.dat" u 1:(1e3*$7) t "Step 1" ls 1, \ - folder2."cpuTime.dat" u 1:(1e3*$7) t "Step 2" ls 2, \ - folder3."cpuTime.dat" u 1:(1e3*$7) t "Step 3" ls 3, \ - folder4."cpuTime.dat" u 1:(1e3*$7) t "Step 4" ls 4 - \ No newline at end of file +plot folder1."cpuTime.dat" u 1:(1e3*$7) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e3*$7) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e3*$7) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e3*$7) t name4 ls 4 + +set title "Num. particles" +plot folder1."cpuTime.dat" u 1:($2) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:($2) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:($2) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:($2) t name4 ls 4 + +unset multiplot + + +set terminal qt enhanced 2 persist size 1600, 1000 font "Times ,10" + +set style line 1 pt 4 lc rgb "#B50427" #Squares red +set style line 2 pt 6 lc rgb "#3B4CC1" #Circles blue +set style line 3 pt 1 lc rgb "#2CA02C" #Crosses green +set style line 4 pt 2 lc rgb "#FE7F0E" #Exes orange +set style line 5 pt 8 lc rgb "#D6696B" #Triangles light red +set style line 10 lt 1 lw 2.0 lc rgb "black" #Black solid line + +set key box opaque +set pointsize 1.5 + +set ylabel "Time per particle (micros)" +set xrange [0:500] +set xlabel "Iteration" + + +set multiplot layout 2,3 +set title "Total" +plot folder1."cpuTime.dat" u 1:(1e6*$3/$2) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e6*$3/$2) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e6*$3/$2) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e6*$3/$2) t name4 ls 4 + +set title "Push" +plot folder1."cpuTime.dat" u 1:(1e6*$4/$2) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e6*$4/$2) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e6*$4/$2) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e6*$4/$2) t name4 ls 4 + +set title "Reset" +plot folder1."cpuTime.dat" u 1:(1e6*$5/$2) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e6*$5/$2) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e6*$5/$2) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e6*$5/$2) t name4 ls 4 + +set title "Collisions" +plot folder1."cpuTime.dat" u 1:(1e6*$6/$2) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e6*$6/$2) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e6*$6/$2) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e6*$6/$2) t name4 ls 4 + +set title "Weighting" +plot folder1."cpuTime.dat" u 1:(1e6*$7/$2) t name1 ls 1, \ + folder2."cpuTime.dat" u 1:(1e6*$7/$2) t name2 ls 2, \ + folder3."cpuTime.dat" u 1:(1e6*$7/$2) t name3 ls 3, \ + folder4."cpuTime.dat" u 1:(1e6*$7/$2) t name4 ls 4 + +unset multiplot diff --git a/runs/cylFlow/plotCPUTime.gp b/runs/cylFlow/plotCPUTime.gp index 62c9314..58d8ba4 100644 --- a/runs/cylFlow/plotCPUTime.gp +++ b/runs/cylFlow/plotCPUTime.gp @@ -1,5 +1,5 @@ reset -set terminal qt enhanced 1 size 1800, 400 font "Times ,10" +set terminal qt enhanced 1 persist size 1800, 400 font "Times ,10" set style line 1 pt 4 lc rgb "#B50427" #Squares red set style line 2 pt 6 lc rgb "#3B4CC1" #Circles blue @@ -46,4 +46,4 @@ set xlabel "Iteration" plot "cpuTime.dat" u 1:(1e6*$4/$2) t "Push" ls 1, \ "" u 1:(1e6*$5/$2) t "Reset" ls 2, \ "" u 1:(1e6*$6/$2) t "Collisions" ls 3, \ - "" u 1:(1e6*$7/$2) t "Weighting" ls 4 \ No newline at end of file + "" u 1:(1e6*$7/$2) t "Weighting" ls 4 diff --git a/src/DSMC_Neutrals.f95 b/src/DSMC_Neutrals.f95 index 7d6ab1d..310f393 100644 --- a/src/DSMC_Neutrals.f95 +++ b/src/DSMC_Neutrals.f95 @@ -43,9 +43,9 @@ PROGRAM DSMC_Neutrals !$OMP PARALLEL PRIVATE(t) DEFAULT(SHARED) DO t = 1, tmax - tStep = omp_get_wtime() !Insert new particles !$OMP SINGLE + tStep = omp_get_wtime() tPush = omp_get_wtime() DO i=1, nInject CALL inject(i)%addParticles() @@ -55,6 +55,7 @@ PROGRAM DSMC_Neutrals !$OMP DO DO n=1, n_part_old CALL push(part_old(n)) + CALL mesh%vols(part_old(n)%e_p)%obj%findCell(part_old(n)) END DO !$OMP END DO @@ -62,10 +63,12 @@ PROGRAM DSMC_Neutrals !$OMP SINGLE tPush = omp_get_wtime() - tPush tReset = omp_get_wtime() + !$OMP END SINGLE !Reset particles - CALL resetParticles() + CALL resetParticles(part_inj, part_old) + !$OMP SINGLE tReset = omp_get_wtime() - tReset tColl = omp_get_wtime() !$OMP END SINGLE @@ -73,7 +76,6 @@ PROGRAM DSMC_Neutrals !Collisions !$OMP DO DO e=1, mesh%numVols - mesh%vols(e)%obj%nColl = 0 !Reset number of collisions CALL mesh%vols(e)%obj%collision() END DO !$OMP END DO @@ -85,7 +87,7 @@ PROGRAM DSMC_Neutrals tWeight = omp_get_wtime() !$OMP END SINGLE - CALL scatterGrid() + CALL scatterGrid(mesh, part_old) !$OMP SINGLE tWeight = omp_get_wtime() - tWeight diff --git a/src/modules/makefile b/src/modules/makefile index a7df1b5..d8ee4ae 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -16,7 +16,7 @@ moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver. moduleInject.o: moduleSpecies.o moduleSolver.o moduleMesh.o moduleMeshCyl.o moduleInject.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ -moduleList.o: moduleErrors.o moduleList.f95 +moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ moduleMesh.o: moduleOutput.o moduleList.o moduleSpecies.o moduleMesh.f95 @@ -37,7 +37,7 @@ moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95 moduleRefParam.o: moduleConstParam.o moduleRefParam.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ -moduleSpecies.o: moduleCaseParam.o moduleList.o moduleSpecies.f95 +moduleSpecies.o: moduleCaseParam.o moduleSpecies.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95, $@) -o $(OBJDIR)/$@ moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95 diff --git a/src/modules/moduleCollisions.f95 b/src/modules/moduleCollisions.f95 index f59084b..7367b29 100644 --- a/src/modules/moduleCollisions.f95 +++ b/src/modules/moduleCollisions.f95 @@ -17,13 +17,13 @@ MODULE moduleCollisions END SUBROUTINE - SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrel, part_i, part_j) + SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j) USE moduleSpecies IMPORT:: collisionBinary CLASS(collisionBinary), INTENT(in):: self REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(out):: sigmaVrel + REAL(8), INTENT(inout):: sigmaVrelMaxNew TYPE(particle), INTENT(inout):: part_i, part_j END SUBROUTINE @@ -115,7 +115,8 @@ MODULE moduleCollisions END SUBROUTINE - SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrel, part_i, part_j) + SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, & + part_i, part_j) USE moduleConstParam USE moduleSpecies USE moduleTable @@ -123,8 +124,9 @@ MODULE moduleCollisions CLASS(collisionBinaryElastic), INTENT(in):: self REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(out):: sigmaVrel + REAL(8), INTENT(inout):: sigmaVrelMaxNew TYPE(particle), INTENT(inout):: part_i, part_j + REAL(8):: sigmaVrel REAL(8):: vRel !relative velocity REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities REAL(8):: alpha @@ -132,7 +134,8 @@ MODULE moduleCollisions !v relative (in units of [m][L]^2[s]^-2 vRel = species(1)%obj%m*NORM2(part_i%v - part_j%v)**2 sigmaVrel = self%crossSec%get(vRel)*vRel - IF (sigmaVrel/sigmaVrelMax > RAND()) THEN + sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel + IF (sigmaVrelMaxNew/sigmaVrelMax > RAND()) THEN !Applies the collision v_i = NORM2(part_i%v) v_j = NORM2(part_j%v) @@ -147,6 +150,7 @@ MODULE moduleCollisions END IF + END SUBROUTINE END MODULE moduleCollisions diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95 index 7b6ff28..e0fea6d 100644 --- a/src/modules/moduleInject.f95 +++ b/src/modules/moduleInject.f95 @@ -178,6 +178,8 @@ MODULE moduleInject !Push new particle CALL push(part_inj(n)) + !Assign cell to new particle + CALL mesh%vols(part_inj(n)%e_p)%obj%findCell(part_inj(n)) END DO diff --git a/src/modules/moduleList.f95 b/src/modules/moduleList.f95 index 830688b..07a4013 100644 --- a/src/modules/moduleList.f95 +++ b/src/modules/moduleList.f95 @@ -1,9 +1,11 @@ !Linked list of particles MODULE moduleList + USE moduleSpecies IMPLICIT NONE + TYPE lNode - INTEGER:: n = 0 - TYPE(lNode), POINTER:: next => NULL() + TYPE(particle), POINTER:: part => NULL() + TYPE(lNode), POINTER:: next => NULL() END TYPE lNode TYPE listNode @@ -17,46 +19,48 @@ MODULE moduleList END TYPE listNode CONTAINS + !Finalize node !Adds element to list - SUBROUTINE addToList(this,n) - INTEGER,INTENT(in):: n - CLASS(listNode):: this + SUBROUTINE addToList(this,part) + USE moduleSpecies + CLASS(listNode), INTENT(inout):: this + TYPE(particle),INTENT(in), TARGET:: part TYPE(lNode),POINTER:: temp ALLOCATE(temp) - temp%n = n + temp%part => part NULLIFY(temp%next) + this%amount = this%amount + 1 IF (.NOT. ASSOCIATED(this%head)) THEN !First element this%head => temp this%tail => temp - this%amount = 1 ELSE !Append element this%tail%next => temp this%tail => temp - this%amount = this%amount + 1 END IF END SUBROUTINE addToList - FUNCTION getFromList(self, iObj) RESULT(n) + FUNCTION getFromList(self, iObj) RESULT(partTemp) use moduleErrors IMPLICIT NONE CLASS(listNode):: self INTEGER:: iObj - INTEGER:: n + TYPE(particle), POINTER:: partTemp INTEGER:: i - TYPE(lNode), POINTER:: tempNode + TYPE(lNode), POINTER:: nodeTemp + IF (iObj > self%amount) CALL criticalError('Accessing to element list outisde range','getFromList') - tempNode => self%head + nodeTemp => self%head DO i = 1, iObj - 1 - tempNode => tempNode%next + nodeTemp => nodeTemp%next END DO - n = tempNode%n + partTemp => nodeTemp%part END FUNCTION getFromList diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95 index 9d3d3a6..10ea618 100644 --- a/src/modules/moduleMesh.f95 +++ b/src/modules/moduleMesh.f95 @@ -93,10 +93,11 @@ MODULE moduleMesh !Number of collisions per volume INTEGER:: nColl = 0 CONTAINS - PROCEDURE(initVol_interface), DEFERRED, PASS:: init - PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter - PROCEDURE(collision_interface), DEFERRED, PASS:: collision - PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell + PROCEDURE(initVol_interface), DEFERRED, PASS:: init + PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter + PROCEDURE(collision_interface), DEFERRED, PASS:: collision + PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell + PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput END TYPE meshVol @@ -128,12 +129,18 @@ MODULE moduleMesh USE moduleSpecies IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self + CLASS(meshVol), INTENT(inout):: self CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell CLASS(particle), INTENT(inout):: part END SUBROUTINE findCell_interface + SUBROUTINE resetOutput_interface(self) + IMPORT:: meshVol + CLASS(meshVol), INTENT(inout):: self + + END SUBROUTINE resetOutput_interface + END INTERFACE !Containers for volumes in the mesh diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95 index 37fcc05..61b9fb6 100644 --- a/src/modules/moduleMeshCyl.f95 +++ b/src/modules/moduleMeshCyl.f95 @@ -62,17 +62,18 @@ MODULE moduleMeshCyl REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolQuadCyl - PROCEDURE, PASS:: locKe => localKeQuad - PROCEDURE, PASS:: detJac => detJQuad - PROCEDURE, PASS:: invJ => invJQuad - PROCEDURE, PASS:: area => areaQuad - PROCEDURE, NOPASS:: weight => weightQuad - PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: dVal => dValQuad - PROCEDURE, PASS:: scatter => scatterQuad - PROCEDURE, PASS:: phy2log => phy2logQuad - PROCEDURE, PASS:: findCell => findCellCylQuad + PROCEDURE, PASS:: init => initVolQuadCyl + PROCEDURE, PASS:: locKe => localKeQuad + PROCEDURE, PASS:: detJac => detJQuad + PROCEDURE, PASS:: invJ => invJQuad + PROCEDURE, PASS:: area => areaQuad + PROCEDURE, NOPASS:: weight => weightQuad + PROCEDURE, NOPASS:: inside => insideQuad + PROCEDURE, PASS:: dVal => dValQuad + PROCEDURE, PASS:: scatter => scatterQuad + PROCEDURE, PASS:: phy2log => phy2logQuad + PROCEDURE, PASS:: findCell => findCellCylQuad + PROCEDURE, PASS:: resetOutput => resetOutputQuad END TYPE meshVolCylQuad @@ -430,7 +431,7 @@ MODULE moduleMeshCyl USE moduleSpecies IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVolCylQuad), INTENT(inout):: self CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell CLASS(particle), INTENT(inout):: part REAL(8):: xLog(1:2) @@ -482,6 +483,36 @@ MODULE moduleMeshCyl END SUBROUTINE findCellCylQuad + PURE SUBROUTINE resetOutputQuad(self) + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(inout):: self + INTEGER:: k + + DO k = 1, nSpecies + self%n1%output(k)%den = 0.D0 + self%n1%output(k)%mom = 0.D0 + self%n1%output(k)%tensorS = 0.D0 + + self%n2%output(k)%den = 0.D0 + self%n2%output(k)%mom = 0.D0 + self%n2%output(k)%tensorS = 0.D0 + + self%n3%output(k)%den = 0.D0 + self%n3%output(k)%mom = 0.D0 + self%n3%output(k)%tensorS = 0.D0 + + self%n4%output(k)%den = 0.D0 + self%n4%output(k)%mom = 0.D0 + self%n4%output(k)%tensorS = 0.D0 + + END DO + + END SUBROUTINE resetOutputQuad + + SUBROUTINE collision2DCyl(self) USE moduleRefParam USE moduleConstParam @@ -492,45 +523,41 @@ MODULE moduleMeshCyl CLASS(meshVolCyl), INTENT(inout):: self INTEGER:: Npart !Number of particles inside the cell - INTEGER:: Ncoll !Maximum number of collisions REAL(8):: Fn !Specific weight REAL(8):: Pmax !Maximum probability of collision - INTEGER:: i,j !random particles index INTEGER:: rnd !random index TYPE(particle), POINTER:: part_i, part_j INTEGER:: n !collision INTEGER:: ij, k - REAL(8):: sigmaVrel REAL(8):: sigmaVrelMaxNew - Fn = species(1)%obj%weight!Check how to do this for multiple species + Fn = species(1)%obj%weight!TODO: Check how to do this for multiple species Npart = self%listPart_in%amount Pmax = Fn*self%sigmaVrelMax*tau/self%volume - Ncoll = INT(REAL(Npart*(Npart-1))*Pmax*0.5D0) - self%nColl = Ncoll + self%nColl = INT(REAL(Npart*(Npart-1))*Pmax*0.5D0) - DO n = 1, NColl + DO n = 1, self%NColl !Select random numbers rnd = 1 + FLOOR(Npart*RAND()) - i = self%listPart_in%get(rnd) - part_i => part_old(i) + part_i => self%listPart_in%get(rnd) rnd = 1 + FLOOR(Npart*RAND()) - j = self%listPart_in%get(rnd) - part_j => part_old(j) + part_j => self%listPart_in%get(rnd) ij = interactionIndex(part_i%pt, part_j%pt) sigmaVrelMaxNew = 0.D0 DO k = 1, interactionMatrix(ij)%amount - CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrel, part_i, part_j) - sigmaVrelMaxNew = sigmaVrel + SigmaVrelMaxNew + CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) END DO !Update maximum cross section times vrelative per each collision IF (sigmaVrelMaxNew > self%sigmaVrelMax) self%sigmaVrelMax = sigmaVrelMaxNew - END DO + !Reset output in nodes + CALL self%resetOutput() + + !Erase the list of particles inside the cell CALL self%listPart_in%erase() END SUBROUTINE collision2DCyl diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95 index 48f3e7e..bf23dc5 100644 --- a/src/modules/moduleSolver.f95 +++ b/src/modules/moduleSolver.f95 @@ -2,38 +2,25 @@ MODULE moduleSolver CONTAINS - SUBROUTINE scatterGrid() + SUBROUTINE scatterGrid(meshIn, partArray) USE moduleSpecies - USE moduleRefParam USE moduleMesh - USE moduleOutput - INTEGER:: n, e, k - !Cleans previous output - !$OMP DO PRIVATE(k) - DO e = 1, mesh%numNodes - DO k= 1, nSpecies - mesh%nodes(e)%obj%output(k)%den = 0.D0 - mesh%nodes(e)%obj%output(k)%mom = 0.D0 - mesh%nodes(e)%obj%output(k)%tensorS = 0.D0 - - END DO - - END DO - !$OMP END DO + CLASS(meshGeneric), INTENT(in):: meshIn + TYPE(particle), INTENT(in):: partArray(:) + INTEGER:: n !Loops over the particles to scatter them !$OMP DO DO n=1, n_part_old - CALL mesh%vols(part_old(n)%e_p)%obj%scatter(part_old(n)) + CALL meshIn%vols(partArray(n)%e_p)%obj%scatter(partArray(n)) END DO !$OMP END DO - END SUBROUTINE scatterGrid - SUBROUTINE push(part) + PURE SUBROUTINE push(part) USE moduleSpecies USE moduleMesh IMPLICIT NONE @@ -68,44 +55,66 @@ MODULE moduleSolver part_temp%e_p = part%e_p !Assign cell to particle part=part_temp - CALL mesh%vols(part%e_p)%obj%findCell(part) END SUBROUTINE push - SUBROUTINE resetParticles() + SUBROUTINE resetParticles(partInj, partOld) USE moduleSpecies USE moduleList USE moduleMesh IMPLICIT NONE + TYPE(particle), INTENT(in), ALLOCATABLE:: partInj(:) + TYPE(particle), INTENT(inout), ALLOCATABLE:: partOld(:) INTEGER:: nn, n + INTEGER, SAVE:: n_inj_in, n_old_in TYPE(particle), ALLOCATABLE:: partTemp(:) - IF (n_part_old > 0) THEN - n_part_new = COUNT(part_old%n_in) + COUNT(part_inj%n_in) - ELSE - n_part_new = COUNT(part_inj%n_in) - END IF + !$OMP SECTIONS + !$OMP SECTION + n_inj_in = 0 + IF (ALLOCATED(partInj)) THEN + n_inj_in = COUNT(partInj%n_in) + END IF + !$OMP SECTION + n_old_in = 0 + IF (ALLOCATED(partOld)) THEN + n_old_in = COUNT(partOld%n_in) + + END IF + !$OMP END SECTIONS + + !$OMP BARRIER + + !$OMP SINGLE CALL MOVE_ALLOC(part_old, partTemp) - ALLOCATE(part_old(1:n_part_new)) + n_part_new = n_inj_in + n_old_in + ALLOCATE(partOld(1:n_part_new)) nn = 0 DO n = 1, nPartInj IF (part_inj(n)%n_in) THEN nn = nn + 1 part_old(nn) = part_inj(n) - CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(nn) + CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(part_old(nn)) + END IF + END DO + DO n = 1, n_part_old IF (partTemp(n)%n_in) THEN nn = nn + 1 - part_old(nn) = partTemp(n) - CALL mesh%vols(part_old(nn)%e_p)%obj%listPart_in%add(nn) + partOld(nn) = partTemp(n) + CALL mesh%vols(partOld(nn)%e_p)%obj%listPart_in%add(partOld(nn)) + END IF + END DO + n_part_old = n_part_new + !$OMP END SINGLE END SUBROUTINE resetParticles diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95 index 113947e..6ca4b87 100644 --- a/src/modules/moduleSpecies.f95 +++ b/src/modules/moduleSpecies.f95 @@ -1,7 +1,6 @@ !Contains the information about species (particles) MODULE moduleSpecies USE moduleCaseParam - USE moduleList IMPLICIT NONE TYPE, ABSTRACT:: speciesGeneric