From 7fd4de24af91797aed57b889199416e5df84c58b Mon Sep 17 00:00:00 2001
From: holger <holger-grosshans@gmx.de>
Date: Mon, 20 Apr 2020 12:06:50 +0200
Subject: [PATCH] 	modified:   run.sh 	modified:  
 src/electrostatics.f90 	modified:   src/makefile 	modified:  
 src/particles.f90 	modified:   src/particlesTransport.f90

---
 run.sh                     |  2 +-
 src/electrostatics.f90     |  2 +-
 src/makefile               |  4 ++--
 src/particles.f90          | 10 +++++++---
 src/particlesTransport.f90 | 34 +++++++++++++---------------------
 5 files changed, 24 insertions(+), 28 deletions(-)

diff --git a/run.sh b/run.sh
index d24482a..ec03eaf 100644
--- a/run.sh
+++ b/run.sh
@@ -1,2 +1,2 @@
 mkdir -p results restart output
-mpiexec -np 2 src/pafiX
+mpiexec -np 1 src/pafiX
diff --git a/src/electrostatics.f90 b/src/electrostatics.f90
index 746ebd5..67c5e8d 100644
--- a/src/electrostatics.f90
+++ b/src/electrostatics.f90
@@ -215,7 +215,7 @@
       use var
       real(kind=pr) :: dis,dis3,disx,disy,disz, &
           Extot,Eytot,Eztot,partmass
-      integer :: n1,n2,ip,jp,lp,m1,m2,n1end
+      integer :: n1,n2,ip,jp,lp,m1,m2,i,j,l
 
 
       select case (elforceScheme)
diff --git a/src/makefile b/src/makefile
index 148b351..765e26b 100644
--- a/src/makefile
+++ b/src/makefile
@@ -15,9 +15,9 @@ INC =
 # gcc:
 CMP = mpifort
 #O3: optimization
-FLAGS   = -O3 -mcmodel=medium -pg
+#FLAGS   = -O3 -mcmodel=medium -pg
 #debugging:
-#FLAGS   = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
+FLAGS   = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
 FFLAGS  = -c $(FLAGS)
 
 
diff --git a/src/particles.f90 b/src/particles.f90
index 6e53b44..eb9017b 100644
--- a/src/particles.f90
+++ b/src/particles.f90
@@ -3,6 +3,7 @@
 !> @brief solve particulate phase
       subroutine solveParticles
       use var
+      integer ip,jp,lp
 
       if (pnd.ne.0.or.npTot.ne.0) then
         if (myid.eq.0) write(*,'(a)',advance='no') 'particle'
@@ -347,11 +348,13 @@
 !> @brief locate particles on Eulerian grid
       subroutine particlesToCells
       use var
+      integer :: i,j,l,n,ip,jp,lp
 
-      npic=0
+      npic(:,:,:)=0
 
       do 1 n=1,np
-        npic(ip(n),jp(n),lp(n))=npic(ip(n),jp(n),lp(n))+1
+        i=ip(n); j=jp(n); l=lp(n)
+        npic(i,j,l)=npic(i,j,l)+1
 1     enddo
 
       if (allocated(nic)) deallocate(nic)
@@ -360,7 +363,8 @@
       nic=0
 
       do 2 n=1,np
-        nic(ip(n),jp(n),lp(n),count(nic(ip(n),jp(n),lp(n),:).ne.0)+1)=n
+        i=ip(n); j=jp(n); l=lp(n)
+        nic(i,j,l,count(nic(i,j,l,:).ne.0)+1)=n
 2     enddo
 
 
diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90
index 36c31b2..7105231 100644
--- a/src/particlesTransport.f90
+++ b/src/particlesTransport.f90
@@ -59,7 +59,7 @@
         cd= 24._pr/Reyp * (1._pr + 1._pr/6._pr*Reyp**(2._pr/3._pr))
       endif
       drag= 0.375_pr*rhof*urel*cd/(rhop*radp(n))
-!      drag=0._pr
+      !drag=0._pr
 
 
       return
@@ -170,6 +170,7 @@
           if ((zp(n).lt.zmin+radp(n)).or.(zp(n).gt.zmax-radp(n))) then
             wp(n)=-wp(n)*restRatio      
             zp(n)=zpold+wp(n)*dt
+            wcollnum(n)=wcollnum(n)+1
             if (qpmax.gt.qp0) call chargeParticleWall(n,3)
           endif
         endif
@@ -356,7 +357,7 @@
                        r3n1,r3n2,rxx,projvtox,srel,projvton, &
                        closestrelvelo,fracdt, &
                        vex,vey,vez,vnx,vny,vnz
-      integer :: n1,n2,numcol,ip,jp,lp,m1,m2
+      integer :: n1,n2,numcol,ip,jp,lp,m1,m2,i,j,l
       
 
       numcol=0
@@ -364,16 +365,13 @@
       do 3 i=1,ii
       do 3 j=1,jj
       do 3 l=1,ll
+! condition I only particles in same cell:
       do 1 m1=2,npic(i,j,l)
       do 2 m2=1,(m1-1)
 
       n1=nic(i,j,l,m1)
       n2=nic(i,j,l,m2)
 
-! condition I to check collision (same or neighbour cell):
-      if ((ip(n1).lt.ip(n2)-1.or.ip(n1).gt.ip(n2)+1).or. &
-          (jp(n1).lt.jp(n2)-1.or.jp(n1).gt.jp(n2)+1).or. &
-          (lp(n1).lt.lp(n2)-1.or.lp(n1).gt.lp(n2)+1)) cycle
 ! setting the 2 particles in the rest frame of particle n2
       relvelo(1)=up(n1)-up(n2)
       relvelo(2)=vp(n1)-vp(n2)
@@ -391,17 +389,17 @@
 ! condition IIb to check collision (no rel velocity, no collision):
       if (abs(srel).lt.1.e-19_pr) cycle
 
-      partdist=dot_product(midlink,midlink)
+      partdist=sqrt(dot_product(midlink,midlink))
       sumrad=radp(n1)+radp(n2)
       
 ! condition III to check collision (contact expected):
-      if ((partdist-(projvtox/srel)**2).gt.sumrad**2) then
+      if ((partdist**2-(projvtox/srel)**2).gt.sumrad**2) then
         cycle
       elseif (partdist.le.sumrad) then ! if overlapping after seeding
         cycle
       endif
 
-      closestrelvelo=sqrt(sumrad**2-(partdist-(projvtox/srel)**2))
+      closestrelvelo=sqrt(sumrad**2-(partdist**2-(projvtox/srel)**2))
       fracdt=(projvtox/srel-closestrelvelo)/srel
 
 ! condition IV to check collision (contact during next timestep):
@@ -510,7 +508,7 @@
         endif
       else
         npseed=int(pnd*dimx*dimy*dimz)
-!       np=2
+ !      npseed=2
       endif
       npseedTot=syncSumI(npseed)
       np=npp+npseed
@@ -545,17 +543,6 @@
         yp(n)=randomy(n-npp)*(dimy-2._pr*prad)+yf(jmin-1)+prad
         zp(n)=randomz(n-npp)*(dimz-2._pr*prad)+zf(lmin-1)+prad
 
-!      if (mod(n,4).eq.0) yp(n)=-1.99e-2_pr
-!      if (mod(n,4).eq.1) yp(n)=1.99e-2_pr
-!        yp(n)=0_pr
-!      if (mod(n,4).eq.2) zp(n)=-1.99e-2_pr
-!      if (mod(n,4).eq.3) zp(n)=1.99e-2_pr
-
-! yp(1)=-0.019_pr;zp(1)=-0.019_pr
-! yp(2)=0.019_pr;zp(2)=0.019_pr
-! vp(1)=-1._pr;vp(2)=1._pr;wp(1)=-1._pr;wp(2)=1._pr
-
-! if(myid.eq.0) print*,xp(n),yp(n),zp(n),up(n),vp(n),wp(n)
  
         wcollnum(n)=0
         ppcollnum(n)=0
@@ -571,6 +558,11 @@
         wp(n)=w(i,j,l)
         if (zp(n).lt.zc(l)) wp(n)=w(i,j,l-1)
       
+
+!      xp(1)=0.02;yp(1)=0.;zp(1)=0.
+!      xp(2)=0.03;yp(2)=0.;zp(2)=0.
+!      up(1)=1.;up(2)=-1;vp(1)=0.;vp(2)=0.;wp(1)=0.;wp(2)=0.
+
         uf(n)=up(n); vf(n)=vp(n); wf(n)=wp(n);
         radp(n)=prad
         partn(n)=1._pr
-- 
GitLab