From 27d4b72a4a2e514abd96d4d8d3cd6f707354c224 Mon Sep 17 00:00:00 2001 From: holger <holger-grosshans@gmx.de> Date: Sun, 9 Aug 2020 10:06:38 +0200 Subject: [PATCH] modified: src/bc.f90 modified: src/electrostatics.f90 modified: src/main.f90 modified: src/parallel.f90 modified: src/particles.f90 modified: src/particlesTransport.f90 modified: src/pre.f90 modified: src/restart.f90 modified: src/var.f90 modified: src/writevtk_particles.f90 --- src/bc.f90 | 7 +- src/electrostatics.f90 | 32 ++++---- src/main.f90 | 4 +- src/parallel.f90 | 20 ++--- src/particles.f90 | 146 +++++++++++++++++-------------------- src/particlesTransport.f90 | 16 +++- src/pre.f90 | 2 +- src/restart.f90 | 7 +- src/var.f90 | 8 +- src/writevtk_particles.f90 | 2 +- 10 files changed, 127 insertions(+), 117 deletions(-) diff --git a/src/bc.f90 b/src/bc.f90 index d323585..a211c92 100644 --- a/src/bc.f90 +++ b/src/bc.f90 @@ -3,7 +3,6 @@ !> @brief boundary condition for velocity and pressure field subroutine bcUVWP use var - real(kind=pr),dimension(dimj*diml) :: random real(kind=pr) :: & ue,uw,ve,vw,we,ww, & uba,ufr,vba,vfr,wba,wfr, & @@ -119,9 +118,9 @@ if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j)) if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l)) do m=1,gc - u(imin-m,j,l)=ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk0 -! v(imin-m,j,l)=0._pr -! w(imin-m,j,l)=0._pr + u(imin-m,j,l)=max(0._pr,ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk0) + v(imin-m,j,l)=0._pr + w(imin-m,j,l)=0._pr enddo enddo; enddo diff --git a/src/electrostatics.f90 b/src/electrostatics.f90 index 39928c3..142a19f 100644 --- a/src/electrostatics.f90 +++ b/src/electrostatics.f90 @@ -7,7 +7,7 @@ integer it - if (pnd.ne.0._pr.and.syncMax(maxval(q_el)).ne.0._pr) then + if (pnd.ne.0._pr.and.syncMax(maxval(abs(q_el))).ne.0._pr) then do it=1,itmax call calcPhi_el(err) @@ -283,7 +283,7 @@ !> @brief compute particle-wall charging subroutine chargeParticleWall(n,direction) use var - real(kind=pr) :: ut2,un2,alpha1,AoAtot,dqp + real(kind=pr) :: ut2,un2,alpha1,AoAtot,dqp,random(2),random_normal integer :: n,i,direction character(70) :: filename @@ -307,16 +307,18 @@ ! alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*un2* & ! ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr) ! AoAtot=alpha1/4._pr - AoAtot=1._pr +! dqp=Qaccfactor*AoAtot*(qpmax-q_el(n)) - dqp=Qaccfactor*AoAtot*(qpmax-q_el(n)) - q_el(n)=min(dqp+q_el(n),qpmax) + call random_number(random) + random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method + dqp=Qaccfactor*(qpmax-q_el(n))*(random_normal+1._pr) ! ...(mu*random_normal+sigma) + q_el(n)=max(min(q_el(n)+dqp,qpmax),qp0) write(filename,'(a,i3.3,a)') 'results/particlesWall_p',myid,'.dat' open(access='append',unit=10,file=filename) - write(10,77) nt,t,sqrt(ut2),sqrt(un2),AoAtot,dqp,q_el(n) -77 format(i7,6(x,es11.2e2)) + write(10,77) nt,t,sqrt(ut2),sqrt(un2),dqp,q_el(n) +77 format(i7,5(x,es11.3e2)) close(10) return @@ -328,7 +330,7 @@ !> @brief compute particle-particle charging subroutine chargeParticleParticle(n1,n2) use var - real(kind=pr) urel2,alpha1,dqp + real(kind=pr) urel2,alpha1,dqp,AoAtot,random(2),random_normal integer :: n1,n2 character(70) :: filename @@ -343,19 +345,19 @@ ! *(1._pr-nyp**2)/Ep)**(0.4_pr) ! dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1)) ! endif +! dqp=Qaccfactor*AoAtot*(q_el(n2)-q_el(n1)) - AoAtot=1._pr - - dqp=Qaccfactor*AoAtot*(q_el(n2)-q_el(n1)) - - q_el(n1)=q_el(n1)+dqp - q_el(n2)=q_el(n2)-dqp + call random_number(random) + random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method + dqp=Qaccfactor*(q_el(n2)-q_el(n1))*(random_normal+1._pr) + q_el(n1)=max(min(q_el(n1)+dqp,qpmax),qp0) + q_el(n2)=max(min(q_el(n2)-dqp,qpmax),qp0) write(filename,'(a,i3.3,a)') 'results/particlesPart_p',myid,'.dat' open(access='append',unit=10,file=filename) write(10,77) nt,t,q_el(n1),q_el(n2),dqp -77 format(i7,5(x,es11.2e2)) +77 format(i7,4(x,es11.3e2)) close(10) diff --git a/src/main.f90 b/src/main.f90 index 45fa4af..91feab5 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -25,11 +25,11 @@ call postProcessing - timecom=syncAv(timecom) + timecom(1)=syncAv(timecom) if (myid.eq.0) write(*,'(3(a,es8.2e2))') & 'time (s) |physical= ',t, & ' |comput. = ',(mpi_wtime()-timebeg), & - ' |MPI com = ',timecom + ' |MPI com = ',timecom(1) 10 enddo diff --git a/src/parallel.f90 b/src/parallel.f90 index 32bf875..8994760 100644 --- a/src/parallel.f90 +++ b/src/parallel.f90 @@ -14,7 +14,7 @@ timenow=real(mpi_wtime(),kind=pr) timebeg=real(mpi_wtime(),kind=pr) timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -89,7 +89,7 @@ 3 enddo timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -164,7 +164,7 @@ 3 enddo timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -226,7 +226,7 @@ endif timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -248,7 +248,7 @@ (myvar,syncMax,1,mpi_pr,mpi_max,mpi_comm_world,mpierr) timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -272,7 +272,7 @@ (myvar,syncSum,1,mpi_pr,mpi_sum,mpi_comm_world,mpierr) timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end function @@ -295,7 +295,7 @@ (myvar,syncSumI,1,mpi_integer,mpi_sum,mpi_comm_world,mpierr) timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -319,7 +319,7 @@ syncAv=myvar/nrprocs timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -351,7 +351,7 @@ endif timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end @@ -384,7 +384,7 @@ endif timeend=real(mpi_wtime(),kind=pr) - timecom=timecom+timeend-timenow + timecom(1)=timecom(1)+timeend-timenow return end diff --git a/src/particles.f90 b/src/particles.f90 index 7b1586d..9463709 100644 --- a/src/particles.f90 +++ b/src/particles.f90 @@ -8,7 +8,6 @@ if (pnd.ne.0.or.npTot.ne.0) then if (myid.eq.0) write(*,'(a)',advance='no') 'particle' if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed - call particlesToCells call fluidVelocity if (elForceScheme.ne.2) call forcesGauss if (elForceScheme.ne.1) call forcesCoulomb @@ -16,7 +15,7 @@ call particlesCollide call particlesTransport !particles from neighbour gc added call chargeDensity - call momentumCoupling !particles from neighbour gc removed + call momentumCoupling !particles from neighbour gc removed if (myid.eq.0) write(*,*) endif @@ -49,6 +48,7 @@ partn(m)=partn(n) wcollnum(m)=wcollnum(n) ppcollnum(m)=ppcollnum(n) + nGlob(m)=nGlob(n) return end @@ -168,7 +168,7 @@ lend=lpa endif -! volume on the Eulerien grid +! volume on the Eulerian grid if (direction.eq.0) volE=(xf(iend)-xf(ibeg-1))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1)) if (direction.eq.1) volE=(xc(iend+1)-xc(ibeg))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1)) if (direction.eq.2) volE=(xf(iend)-xf(ibeg-1))*(yc(jend+1)-yc(jbeg))*(zf(lend)-zf(lbeg-1)) @@ -236,7 +236,7 @@ jbeg=jpa; jend=jpa lbeg=lpa; lend=lpa -! volume on the Eulerien grid +! volume on the Eulerian grid if (direction.eq.0) volE=(xf(ipa)-xf(ipa-1))*(yf(jpa)-yf(jpa-1))*(zf(lpa)-zf(lpa-1)) if (direction.eq.1) volE=(xc(ipa+1)-xc(ipa))*(yf(jpa)-yf(jpa-1))*(zf(lpa)-zf(lpa-1)) if (direction.eq.2) volE=(xf(ipa)-xf(ipa-1))*(yc(jpa+1)-yc(jpa))*(zf(lpa)-zf(lpa-1)) @@ -247,7 +247,6 @@ return end - !#################################################################### !> @author Holger Grosshans !> @brief allocate particle arrays @@ -255,7 +254,7 @@ use var real(kind=pr), allocatable, dimension(:) :: & tup,tvp,twp,tuf,tvf,twf,tuf01,tvf01,twf01,txp,typ,tzp,tradp, & - tq_el,tfx_el,tfy_el,tfz_el,tpartn + tq_el,tfx_el,tfy_el,tfz_el,tpartn,tnGlob integer, allocatable, dimension(:) :: & twcollnum,tppcollnum integer :: syncSumI,nppTot @@ -264,80 +263,69 @@ nppTot=syncSumI(npp) maxnp=npTot*(nrprocs)*3 !if nrprocs=1, one particle can be in domain, gc and sync storage - if (nppTot.ge.1) then - allocate( & - tup(maxnp),tvp(maxnp),twp(maxnp), & - tuf(maxnp),tvf(maxnp),twf(maxnp), & - tuf01(maxnp),tvf01(maxnp),twf01(maxnp), & - txp(maxnp),typ(maxnp),tzp(maxnp), & - tradp(maxnp),tq_el(maxnp), & - tfx_el(maxnp),tfy_el(maxnp),tfz_el(maxnp), & - tpartn(maxnp)) - allocate(twcollnum(maxnp),tppcollnum(maxnp)) - - tup(1:npp)=up(1:npp) - tvp(1:npp)=vp(1:npp) - twp(1:npp)=wp(1:npp) - tuf(1:npp)=uf(1:npp) - tvf(1:npp)=vf(1:npp) - twf(1:npp)=wf(1:npp) - tuf01(1:npp)=uf01(1:npp) - tvf01(1:npp)=vf01(1:npp) - twf01(1:npp)=wf01(1:npp) - txp(1:npp)=xp(1:npp) - typ(1:npp)=yp(1:npp) - tzp(1:npp)=zp(1:npp) - tradp(1:npp)=radp(1:npp) - tq_el(1:npp)=q_el(1:npp) - tfx_el(1:npp)=fx_el(1:npp) - tfy_el(1:npp)=fy_el(1:npp) - tfz_el(1:npp)=fz_el(1:npp) - tpartn(1:npp)=partn(1:npp) - twcollnum(1:npp)=wcollnum(1:npp) - tppcollnum(1:npp)=ppcollnum(1:npp) - endif - - if (allocated(up)) then - deallocate( & - up,vp,wp,uf,vf,wf,uf01,vf01,wf01,xp,yp,zp,radp,q_el, & - fx_el,fy_el,fz_el,partn) - deallocate(wcollnum,ppcollnum) - endif + call allocateParticleArray(nppTot,up) + call allocateParticleArray(nppTot,vp) + call allocateParticleArray(nppTot,wp) + call allocateParticleArray(nppTot,uf) + call allocateParticleArray(nppTot,vf) + call allocateParticleArray(nppTot,wf) + call allocateParticleArray(nppTot,uf01) + call allocateParticleArray(nppTot,vf01) + call allocateParticleArray(nppTot,wf01) + call allocateParticleArray(nppTot,xp) + call allocateParticleArray(nppTot,yp) + call allocateParticleArray(nppTot,zp) + call allocateParticleArray(nppTot,radp) + call allocateParticleArray(nppTot,q_el) + call allocateParticleArray(nppTot,fx_el) + call allocateParticleArray(nppTot,fy_el) + call allocateParticleArray(nppTot,fz_el) + call allocateParticleArray(nppTot,partn) + + call allocateParticleArrayI(nppTot,wcollnum) + call allocateParticleArrayI(nppTot,ppcollnum) + call allocateParticleArrayI(nppTot,nGlob) + + contains + + subroutine allocateParticleArray(nppTot,myvar) + use var + real(kind=pr), allocatable, dimension(:) :: & + myvar,tmyvar + integer :: nppTot + + if (nppTot.ge.1) then + allocate(tmyvar(maxnp)) + tmyvar(1:npp)=myvar(1:npp) + endif + if (allocated(myvar)) deallocate(myvar) + allocate(myvar(maxnp)) + myvar=0._pr + if (nppTot.ge.1) then + myvar(1:npp)=tmyvar(1:npp) + endif + + end + + subroutine allocateParticleArrayI(nppTot,myvar) + use var + integer, allocatable, dimension(:) :: & + myvar,tmyvar + integer :: nppTot + + if (nppTot.ge.1) then + allocate(tmyvar(maxnp)) + tmyvar(1:npp)=myvar(1:npp) + endif + if (allocated(myvar)) deallocate(myvar) + allocate(myvar(maxnp)) + myvar=0 + if (nppTot.ge.1) then + myvar(1:npp)=tmyvar(1:npp) + endif + + end - allocate( & - up(maxnp),vp(maxnp),wp(maxnp), & - uf(maxnp),vf(maxnp),wf(maxnp), & - uf01(maxnp),vf01(maxnp),wf01(maxnp), & - xp(maxnp),yp(maxnp),zp(maxnp), & - radp(maxnp),q_el(maxnp), & - fx_el(maxnp),fy_el(maxnp),fz_el(maxnp), & - partn(maxnp)) - allocate(wcollnum(maxnp),ppcollnum(maxnp)) - - if (nppTot.ge.1) then - up(1:npp)=tup(1:npp) - vp(1:npp)=tvp(1:npp) - wp(1:npp)=twp(1:npp) - uf(1:npp)=tuf(1:npp) - vf(1:npp)=tvf(1:npp) - wf(1:npp)=twf(1:npp) - uf01(1:npp)=tuf01(1:npp) - vf01(1:npp)=tvf01(1:npp) - wf01(1:npp)=twf01(1:npp) - xp(1:npp)=txp(1:npp) - yp(1:npp)=typ(1:npp) - zp(1:npp)=tzp(1:npp) - radp(1:npp)=tradp(1:npp) - q_el(1:npp)=tq_el(1:npp) - fx_el(1:npp)=tfx_el(1:npp) - fy_el(1:npp)=tfy_el(1:npp) - fz_el(1:npp)=tfz_el(1:npp) - partn(1:npp)=tpartn(1:npp) - wcollnum(1:npp)=twcollnum(1:npp) - ppcollnum(1:npp)=tppcollnum(1:npp) - endif - - return end diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90 index 6d27428..ee4045f 100644 --- a/src/particlesTransport.f90 +++ b/src/particlesTransport.f90 @@ -204,9 +204,13 @@ call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,uf) call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,vf) call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,wf) + call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,uf01) + call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,vf01) + call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,wf01) call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,partn) call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,wcollnum) call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,ppcollnum) + call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,nGlob) if ((bcx.eq.'p').and.(myid.eq.nrprocs-1)) then @@ -237,6 +241,8 @@ np=m endif + call particlesToCells + return end @@ -262,7 +268,10 @@ do 1 n=1,np - if (ip(n).lt.(imin-1).or.ip(n).gt.(imax+1)) cycle + if ((ip(n).lt.(imin-1)).or.(ip(n).gt.(imax+1)).or. & + (jp(n).lt.(jmin-1)).or.(jp(n).gt.(jmax+1)).or. & + (lp(n).lt.(lmin-1)).or.(lp(n).gt.(lmax+1))) cycle + ! only drag forces f_d(1)= drag(n)*(uf(n)-up(n)) @@ -502,6 +511,8 @@ if (myid.eq.0) write(*,'(a,i9)',advance='no') & ' |seed =',npseedTot + if (npseed.eq.0) return + write(filename,'(a,i3.3,a)') 'results/particlesSeed_p',myid,'.dat' open(access='append',unit=10,file=filename) @@ -528,6 +539,8 @@ 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 + nseedLoc=nseedLoc+1 + nGlob(n)=nseedLoc*1000+myid wcollnum(n)=0 ppcollnum(n)=0 @@ -557,6 +570,7 @@ 1 enddo + call particlesToCells close(10) diff --git a/src/pre.f90 b/src/pre.f90 index d95bb5e..99c0b01 100644 --- a/src/pre.f90 +++ b/src/pre.f90 @@ -193,7 +193,7 @@ ! particles allocate(npic(ii,jj,ll)) - np=0;npp=0 + np=0; npp=0; nseedLoc=0; call allocateParticleArrays ! init diff --git a/src/restart.f90 b/src/restart.f90 index ad3142f..e6a8cbf 100644 --- a/src/restart.f90 +++ b/src/restart.f90 @@ -69,7 +69,12 @@ ntstart=ntstart+1 ntend=ntend+1 - + + do n=1,np + nseedLoc=nseedLoc+1 + nGlob(n)=nseedLoc*1000+myid + enddo + call particlesToCells return diff --git a/src/var.f90 b/src/var.f90 index ee6dd1c..20090b9 100644 --- a/src/var.f90 +++ b/src/var.f90 @@ -118,7 +118,8 @@ np, & !< local number of particles npTot, & !< total number of particles npp, & !< local number of particles before some operation - maxnp !< maximum possible local number of particles + maxnp, & !< maximum possible local number of particles + nseedLoc !< local counter for seeded particles real(kind=pr), allocatable, dimension(:) :: & up,vp,wp, & !< particle velocity (m/s) uf,vf,wf, & !< average fluid velocity around a particle (m/s) @@ -130,7 +131,8 @@ partn !< number of particles represented by parcel integer, allocatable, dimension(:) :: & wcollnum, & !< number of particle-wall collisions - ppcollnum !< number of particle-particle collisions + ppcollnum, & !< number of particle-particle collisions + nGlob !< global particle id integer, allocatable, dimension(:,:,:) :: & npic !< number of particles in a Eulerian cell integer, allocatable, dimension(:,:,:,:) :: & @@ -138,7 +140,7 @@ ! parallel - real(kind=pr) :: timebeg,timecom,timenow,timeend + real(kind=pr) :: timebeg,timecom(10),timenow,timeend integer :: & mpi_pr, & myid, & !< id of current processor diff --git a/src/writevtk_particles.f90 b/src/writevtk_particles.f90 index b61b9d1..b519f06 100644 --- a/src/writevtk_particles.f90 +++ b/src/writevtk_particles.f90 @@ -109,7 +109,7 @@ write(10,'(a)') 'SCALARS np FLOAT 1' write(10,'(a)') 'LOOKUP_TABLE default ' - write(10,fmt=rowfm2) (n,n=1,np) + write(10,fmt=rowfm2) (nGlob(n),n=1,np) write(10,'(a)') write(10,'(a)') 'VECTORS UPVPWP FLOAT' -- GitLab