From 504c97b7d8ad88ae5338146cb00744e6d8187f29 Mon Sep 17 00:00:00 2001 From: holger <holger-grosshans@gmx.de> Date: Sat, 2 May 2020 14:51:56 +0200 Subject: [PATCH] fixed bug with grid mapping derivatives close to the wall --- run.sh | 2 +- src/bc.f90 | 26 +++--- src/electrostatics.f90 | 39 ++++----- src/makefile | 4 +- src/mom5.f90 | 1 + src/momentum.f90 | 168 ++++++++++++++++++------------------- src/parallel.f90 | 151 +++++---------------------------- src/particles.f90 | 6 +- src/particlesTransport.f90 | 78 +++++++---------- src/post.f90 | 2 +- src/pre.f90 | 40 ++++----- src/pressure.f90 | 81 +++++++++--------- src/restart.f90 | 4 +- src/var.f90 | 4 +- src/writedat_fluid.f90 | 6 +- src/writedat_particles.f90 | 2 +- src/writevtk_fluid_xyz.f90 | 2 +- src/writevtk_fluid_xz.f90 | 2 +- src/writevtk_fluid_yz.f90 | 2 +- src/writevtk_grid.f90 | 2 +- src/writevtk_particles.f90 | 2 +- 21 files changed, 237 insertions(+), 387 deletions(-) diff --git a/run.sh b/run.sh index d24482a..aa8cf58 100644 --- a/run.sh +++ b/run.sh @@ -1,2 +1,2 @@ mkdir -p results restart output -mpiexec -np 2 src/pafiX +mpiexec -np 2 --oversubscribe src/pafiX diff --git a/src/bc.f90 b/src/bc.f90 index bf2432b..091ea2e 100644 --- a/src/bc.f90 +++ b/src/bc.f90 @@ -1,7 +1,7 @@ !#################################################################### !> @author Holger Grosshans -!> @brief boundary condition for velocity field - subroutine bcUVW +!> @brief boundary condition for velocity and pressure field + subroutine bcUVWP use var real(kind=pr),dimension(dimj*diml) :: random real(kind=pr) :: & @@ -44,17 +44,17 @@ call random_number(random) alpha=.5_pr mm=0 - do 7 j=jmin,jmax - do 7 l=lmin,lmax - mm=mm+1 - if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l))) - 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 7 m=1,gc - u(imin-m,j,l)=1._pr*ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk - v(imin-m,:,:)=0._pr - w(imin-m,:,:)=0._pr -7 enddo + do j=jmin,jmax; do l=lmin,lmax + mm=mm+1 + if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l))) + 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)=1._pr*ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk + v(imin-m,:,:)=0._pr + w(imin-m,:,:)=0._pr + enddo + enddo; enddo endif if (myid.eq.nrprocs-1) then do 8 m=1,gc diff --git a/src/electrostatics.f90 b/src/electrostatics.f90 index 67c5e8d..47a3fb3 100644 --- a/src/electrostatics.f90 +++ b/src/electrostatics.f90 @@ -44,12 +44,9 @@ err=0._pr phi_el1=0._pr + do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax - do 1 i=imin,imax,1 - do 1 j=jmin,jmax,1 - do 1 l=lmin,lmax,1 - - if (celltype(i,j,l).ne.active) cycle + if (celltype(i,j,l).ne.active) cycle ! solve the poisson equation: diag= (2._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) & @@ -70,7 +67,7 @@ dphi_el=phi_el1(i,j,l)-phi_el(i,j,l) err=err+dphi_el*dphi_el -1 enddo + enddo; enddo; enddo phi_el=phi_el1 @@ -97,11 +94,9 @@ err=0._pr - do 1 i=imin,imax,1 - do 1 j=jmin,jmax,1 - do 1 l=lmin,lmax,1 + do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax - if (celltype(i,j,l).ne.active) cycle + if (celltype(i,j,l).ne.active) cycle ! calculate electrostatic field temp= Ex_el(i,j,l) @@ -118,7 +113,7 @@ err=err+dexel*dexel+deyel*deyel+dezel*dezel -1 enddo + enddo; enddo; enddo call sync(Ex_el); call sync(Ey_el); call sync(Ez_el) call bcE_el @@ -149,14 +144,12 @@ call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0) - do 5 i=ibeg,iend - do 5 j=jbeg,jend - do 5 l=lbeg,lend + do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend iw=i+1-ibeg jw=j+1-jbeg lw=l+1-lbeg rho_el(i,j,l) = rho_el(i,j,l) + q_el(n)*partn(n)*weight(iw,jw,lw)/volE -5 enddo + enddo; enddo; enddo 1 enddo @@ -186,16 +179,14 @@ call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0) - do 5 i=ibeg,iend - do 5 j=jbeg,jend - do 5 l=lbeg,lend + do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend iw=i+1-ibeg jw=j+1-jbeg lw=l+1-lbeg Extot=Extot+Ex_el(i,j,l)*weight(iw,jw,lw) Eytot=Eytot+Ey_el(i,j,l)*weight(iw,jw,lw) Eztot=Eztot+Ez_el(i,j,l)*weight(iw,jw,lw) -5 enddo + enddo; enddo; enddo partmass=4._pr/3._pr*pi*rhop*partn(n)*radp(n)**3 fx_el(n)= q_el(n)*Extot/partmass @@ -232,9 +223,7 @@ 1 enddo case (3) - do 3 i=1,ii - do 3 j=1,jj - do 3 l=1,ll + do i=1,ii; do j=1,jj; do l=1,ll do 4 m1=2,npic(i,j,l) do 5 m2=1,(m1-1) @@ -243,7 +232,7 @@ call forcesCoulombN1N2(n1,n2) 5 enddo 4 enddo -3 enddo + enddo; enddo; enddo end select @@ -296,7 +285,7 @@ use var real(kind=pr) :: ut2,un2,alpha1,AoAtot,dqp integer :: n,i,direction - character*70 :: filename + character(70) :: filename ! pipe: @@ -336,7 +325,7 @@ use var real(kind=pr) uax2,urad2,alpha1,dqp integer :: n1,n2 - character*70 :: filename + character(70) :: filename uax2=(up(n1)-up(n2))**2 urad2=(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2 diff --git a/src/makefile b/src/makefile index 7537dc1..c2b1f4e 100644 --- a/src/makefile +++ b/src/makefile @@ -15,9 +15,9 @@ INC = # gcc: CMP = mpifort #O3: optimization -#FLAGS = -O3 -mcmodel=medium +FLAGS = -O3 -mcmodel=medium #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/mom5.f90 b/src/mom5.f90 index d7ea08e..c2623d0 100644 --- a/src/mom5.f90 +++ b/src/mom5.f90 @@ -42,6 +42,7 @@ call weno(fluxplu,u(i,j,l+3),u(i,j,l+2),u(i,j,l+1),u(i,j,l),u(i,j,l-1)) call weno(fluxmin,u(i,j,l+2),u(i,j,l+1),u(i,j,l),u(i,j,l-1),u(i,j,l-2)) endif +! print*,myid,l,dzcdl(l) wuz= wa*(fluxplu-fluxmin)/dzcdl(l) ! pressure gradient (4th order) diff --git a/src/momentum.f90 b/src/momentum.f90 index 17e394a..5f5366f 100644 --- a/src/momentum.f90 +++ b/src/momentum.f90 @@ -17,47 +17,45 @@ dum=0._pr; dvm=0._pr; dwm=0._pr - do 1 i=imin,imax - do 1 j=jmin,jmax - do 1 l=lmin,lmax - - if (celltype(i,j,l).ne.active) cycle - - dua=0._pr; dva=0._pr; dwa=0._pr - - if (celltype(i+1,j,l).ne.wall) then - call momx1(tu,qu,i,j,l) - du=(defect_u(i,j,l)+tu)/qu*urfu - u1(i,j,l)=u(i,j,l)+du - dua=abs(du) - if (dum.lt.dua) dum=dua - endif - - if (celltype(i,j+1,l).ne.wall) then - call momy1(tv,qv,i,j,l) - dv=(defect_v(i,j,l)+tv)/qv*urfv - v1(i,j,l)=v(i,j,l)+dv - dva=abs(dv) - if (dvm.lt.dva) dvm=dva - endif - - if (celltype(i,j,l+1).ne.wall) then - call momz1(tw,qw,i,j,l) - dw=(defect_w(i,j,l)+tw)/qw*urfw - w1(i,j,l)=w(i,j,l)+dw - dwa=abs(dw) - if (dwm.lt.dwa) dwm=dwa - endif - - -12 err=err+dua*dua+dva*dva+dwa*dwa - -1 enddo + do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax + + if (celltype(i,j,l).ne.active) cycle + + dua=0._pr; dva=0._pr; dwa=0._pr + + if (celltype(i+1,j,l).ne.wall) then + call momx1(tu,qu,i,j,l) + du=(defect_u(i,j,l)+tu)/qu*urfu + u1(i,j,l)=u(i,j,l)+du + dua=abs(du) + if (dum.lt.dua) dum=dua + endif + + if (celltype(i,j+1,l).ne.wall) then + call momy1(tv,qv,i,j,l) + dv=(defect_v(i,j,l)+tv)/qv*urfv + v1(i,j,l)=v(i,j,l)+dv + dva=abs(dv) + if (dvm.lt.dva) dvm=dva + endif + + if (celltype(i,j,l+1).ne.wall) then + call momz1(tw,qw,i,j,l) + dw=(defect_w(i,j,l)+tw)/qw*urfw + w1(i,j,l)=w(i,j,l)+dw + dwa=abs(dw) + if (dwm.lt.dwa) dwm=dwa + endif + + + err=err+dua*dua+dva*dva+dwa*dwa + + enddo; enddo; enddo u=u1; v=v1; w=w1 call sync(u); call sync(v); call sync(w) - call bcUVW + call bcUVWP err=(syncSum(err)/dimgptot/3._pr)**(0.5_pr) ! write(*,'(x,a,es9.2e2,a,3(es9.2e2))') & @@ -100,56 +98,54 @@ ! v(imin:imax, jmin:jmax-1, lmin:lmax) ! w(imin:imax, jmin:jmax, lmin:lmax-1) - do 1 i=imin,imax - do 1 j=jmin,jmax - do 1 l=lmin,lmax - -! if (celltype(i,j,l).ne.active) cycle - - if ((celltype(i+1,j,l).ne.wall).and.(celltype(i-1,j,l).ne.wall).and. & - (celltype(i,j+1,l).ne.wall).and.(celltype(i,j-1,l).ne.wall).and. & - (celltype(i,j,l+1).ne.wall).and.(celltype(i,j,l-1).ne.wall)) then - call mass2(ra,i,j,l) - call mass4(raHO,i,j,l) - defect_c(i,j,l)=raHO-ra - else - defect_c(i,j,l)=0._pr - endif + do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax + +! if (celltype(i,j,l).ne.active) cycle + + if ((celltype(i+1,j,l).ne.wall).and.(celltype(i-1,j,l).ne.wall).and. & + (celltype(i,j+1,l).ne.wall).and.(celltype(i,j-1,l).ne.wall).and. & + (celltype(i,j,l+1).ne.wall).and.(celltype(i,j,l-1).ne.wall)) then + call mass2(ra,i,j,l) + call mass4(raHO,i,j,l) + defect_c(i,j,l)=raHO-ra + else + defect_c(i,j,l)=0._pr + endif - if ((celltype(i+3,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. & - (celltype(i,j+2,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. & - (celltype(i,j,l+2).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then - call momx1(tu,qu,i,j,l) -! call momx3(tuHO,quHO,i,j,l) - call momx5(tuHO,quHO,i,j,l) - defect_u(i,j,l)=tuHO-tu - else - defect_u(i,j,l)=0._pr - endif - - if ((celltype(i+2,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. & - (celltype(i,j+3,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. & - (celltype(i,j,l+2).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then - call momy1(tv,qv,i,j,l) -! call momy3(tvHO,qvHO,i,j,l) - call momy5(tvHO,qvHO,i,j,l) - defect_v(i,j,l)=tvHO-tv - else - defect_v(i,j,l)=0._pr - endif - - if ((celltype(i+2,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. & - (celltype(i,j+2,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. & - (celltype(i,j,l+3).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then - call momz1(tw,qw,i,j,l) -! call momz3(twHO,qwHO,i,j,l) - call momz5(twHO,qwHO,i,j,l) - defect_w(i,j,l)=twHO-tw - else - defect_w(i,j,l)=0._pr - endif - -1 enddo + if ((celltype(i+3,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. & + (celltype(i,j+2,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. & + (celltype(i,j,l+2).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then + call momx1(tu,qu,i,j,l) +! call momx3(tuHO,quHO,i,j,l) + call momx5(tuHO,quHO,i,j,l) + defect_u(i,j,l)=tuHO-tu + else + defect_u(i,j,l)=0._pr + endif + + if ((celltype(i+2,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. & + (celltype(i,j+3,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. & + (celltype(i,j,l+2).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then + call momy1(tv,qv,i,j,l) +! call momy3(tvHO,qvHO,i,j,l) + call momy5(tvHO,qvHO,i,j,l) + defect_v(i,j,l)=tvHO-tv + else + defect_v(i,j,l)=0._pr + endif + + if ((celltype(i+2,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. & + (celltype(i,j+2,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. & + (celltype(i,j,l+3).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then + call momz1(tw,qw,i,j,l) +! call momz3(twHO,qwHO,i,j,l) + call momz5(twHO,qwHO,i,j,l) + defect_w(i,j,l)=twHO-tw + else + defect_w(i,j,l)=0._pr + endif + + enddo; enddo; enddo return diff --git a/src/parallel.f90 b/src/parallel.f90 index fa22848..512fafe 100644 --- a/src/parallel.f90 +++ b/src/parallel.f90 @@ -183,18 +183,18 @@ timenow=real(mpi_wtime(),kind=pr) n=0 - do 1 i = imin,imin+gc-1 - do 1 l = 1,ll; do 1 j = 1,jj + do i= imin,imin+gc-1 + do l= 1,ll; do j= 1,jj n=n+1 - sendleft(n) = myvar(i,j,l) -1 enddo + sendleft(n)= myvar(i,j,l) + enddo; enddo; enddo n=0 - do 2 i = imax,imax-gc+1,-1 - do 2 l = 1,ll; do 2 j = 1,jj + do i= imax,imax-gc+1,-1 + do l= 1,ll; do j= 1,jj n=n+1 - sendright(n) = myvar(i,j,l) -2 enddo + sendright(n)= myvar(i,j,l) + enddo; enddo; enddo call mpi_isend(sendleft,gc*jj*ll,mpi_pr,prev,1,mpi_comm_world,rsleft,mpierr) call mpi_isend(sendright,gc*jj*ll,mpi_pr,next,2,mpi_comm_world,rsright,mpierr) @@ -206,21 +206,21 @@ call mpi_wait(rrleft,mpistatus,mpierr) if (next.ne.mpi_proc_null) then - n=0 - do 3 i = imax+1,ii - do 3 l = 1,ll; do 3 j = 1,jj - n=n+1 - myvar(i,j,l) = recvright(n) -3 enddo + n=0 + do i= imax+1,ii + do l= 1,ll; do j= 1,jj + n=n+1 + myvar(i,j,l)= recvright(n) + enddo; enddo; enddo endif if (prev.ne.mpi_proc_null) then - n=0 - do 4 i = gc,1,-1 - do 4 l = 1,ll; do 4 j = 1,jj - n=n+1 - myvar(i,j,l) = recvleft(n) -4 enddo + n=0 + do i= gc,1,-1 + do l= 1,ll; do j= 1,jj + n=n+1 + myvar(i,j,l)= recvleft(n) + enddo; enddo; enddo endif timeend=real(mpi_wtime(),kind=pr) @@ -245,34 +245,6 @@ call mpi_allreduce& (myvar,syncMax,1,mpi_pr,mpi_max,mpi_comm_world,mpierr) - -! syncMax=myvar -! if (myid.ne.0) then -! call mpi_isend(myvar,1,mpi_pr,0,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -! endif -! -! if (myid.eq.0) then -! do 1 proc=1,nrprocs-1 -! call mpi_irecv(myvar2,1,mpi_pr,proc,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! syncMax=max(myvar,myvar2) -!1 enddo -! endif -! -! if (myid.eq.0) then -! do 2 proc=1,nrprocs-1 -! call mpi_isend(syncMax,1,mpi_pr,proc,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -!2 enddo -! endif -! -! if (myid.ne.0) then -! call mpi_irecv(syncMax,1,mpi_pr,0,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! endif - - timeend=real(mpi_wtime(),kind=pr) timecom=timecom+timeend-timenow @@ -296,34 +268,6 @@ ! mpi_allreduce does not work with intel compiler call mpi_allreduce& (myvar,syncSum,1,mpi_pr,mpi_sum,mpi_comm_world,mpierr) - -! sumvar=myvar -! if (myid.ne.0) then -! call mpi_isend(myvar,1,mpi_pr,0,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -! endif -! -! if (myid.eq.0) then -! do 1 proc=1,nrprocs-1 -! call mpi_irecv(myvar2,1,mpi_pr,proc,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! sumvar=myvar+myvar2 -!1 enddo -! endif -! print*,'x2',myid,sumvar -! -! if (myid.eq.0) then -! do 2 proc=1,nrprocs-1 -! call mpi_isend(sumvar,1,mpi_pr,proc,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -!2 enddo -! endif -! -! if (myid.ne.0) then -! call mpi_irecv(sumvar,1,mpi_pr,0,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! endif -! print*,'x3',myid,myvar timeend=real(mpi_wtime(),kind=pr) timecom=timecom+timeend-timenow @@ -347,33 +291,6 @@ ! mpi_allreduce does not work with intel compiler call mpi_allreduce& (myvar,syncSumI,1,mpi_integer,mpi_sum,mpi_comm_world,mpierr) - - -! syncSumI=myvar -! if (myid.ne.0) then -! call mpi_isend(myvar,1,mpi_integer,0,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -! endif -! -! if (myid.eq.0) then -! do 1 proc=1,nrprocs-1 -! call mpi_irecv(myvar2,1,mpi_integer,proc,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! syncSumI=syncSumI+myvar2 -!1 enddo -! endif -! -! if (myid.eq.0) then -! do 2 proc=1,nrprocs-1 -! call mpi_isend(syncSumI,1,mpi_integer,proc,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -!2 enddo -! endif -! -! if (myid.ne.0) then -! call mpi_irecv(syncSumI,1,mpi_integer,0,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! endif timeend=real(mpi_wtime(),kind=pr) timecom=timecom+timeend-timenow @@ -399,34 +316,6 @@ syncAv=myvar/nrprocs -! syncAv=myvar -! if (myid.ne.0) then -! call mpi_isend(myvar,1,mpi_pr,0,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -! endif -! -! if (myid.eq.0) then -! do 1 proc=1,nrprocs-1 -! call mpi_irecv(myvar2,1,mpi_pr,proc,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! syncAv=syncAv+myvar2 -!1 enddo -! syncAv=syncAv/nrprocs -! endif -! -! if (myid.eq.0) then -! do 2 proc=1,nrprocs-1 -! call mpi_isend(syncAv,1,mpi_pr,proc,0,mpi_comm_world,rs,mpierr) -! call mpi_wait(rs,mpistatus,mpierr) -!2 enddo -! endif -! -! if (myid.ne.0) then -! call mpi_irecv(syncAv,1,mpi_pr,0,0,mpi_comm_world,rr,mpierr) -! call mpi_wait(rr,mpistatus,mpierr) -! endif - - timeend=real(mpi_wtime(),kind=pr) timecom=timecom+timeend-timenow diff --git a/src/particles.f90 b/src/particles.f90 index eb9017b..7b1586d 100644 --- a/src/particles.f90 +++ b/src/particles.f90 @@ -188,9 +188,7 @@ deltaf2=(deltaf**2)/6._pr !compute weights - do 5 i=ibeg,iend - do 5 j=jbeg,jend - do 5 l=lbeg,lend + do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend disx=xp(n)-xc(i) if (direction.eq.1) disx=xp(n)-xf(i) disy=yp(n)-yc(j) @@ -200,7 +198,7 @@ dis2=disx*disx+disy*disy+disz*disz weight(i+1-ibeg,j+1-jbeg,l+1-lbeg)=exp(-dis2/deltaf2) -5 enddo + enddo; enddo; enddo weight=weight/sum(weight) ! sum normalized diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90 index ea3fbd0..7619f49 100644 --- a/src/particlesTransport.f90 +++ b/src/particlesTransport.f90 @@ -12,27 +12,18 @@ uf=0._pr; vf=0._pr; wf=0._pr do 1 n=1,np - if (celltype(ip(n),jp(n),lp(n)).ne.active) cycle - do 4 direction=1,3 - - call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction) - - do 5 i=ibeg,iend - do 5 j=jbeg,jend - do 5 l=lbeg,lend - iw=i+1-ibeg - jw=j+1-jbeg - lw=l+1-lbeg - - if (direction.eq.1) uf(n)=uf(n)+weight(iw,jw,lw)*u(i,j,l) - if (direction.eq.2) vf(n)=vf(n)+weight(iw,jw,lw)*v(i,j,l) - if (direction.eq.3) wf(n)=wf(n)+weight(iw,jw,lw)*w(i,j,l) - -5 enddo + call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction) + do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend + iw=i+1-ibeg + jw=j+1-jbeg + lw=l+1-lbeg + if (direction.eq.1) uf(n)=uf(n)+weight(iw,jw,lw)*u(i,j,l) + if (direction.eq.2) vf(n)=vf(n)+weight(iw,jw,lw)*v(i,j,l) + if (direction.eq.3) wf(n)=wf(n)+weight(iw,jw,lw)*w(i,j,l) + enddo; enddo; enddo 4 enddo - 1 enddo return @@ -113,7 +104,7 @@ np_sendr,np_sendl,np_recvl,np_recvr integer :: n,m,i,j,l integer :: syncSumI - character*70 :: filename + character(70) :: filename logical :: remove(maxnp) remove=.false. @@ -281,31 +272,26 @@ volp= 4._pr/3._pr*pi*partn(n)*radp(n)**3 do 2 direction=1,3 - - call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction) + call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction) ! compute source terms for the parcel - do 5 l=lbeg,lend - do 5 i=ibeg,iend - do 5 j=jbeg,jend - iw=i+1-ibeg - jw=j+1-jbeg - lw=l+1-lbeg - - if (direction.eq.1) then - source= rhop/rhof * volp/volE * f_d(1) - rui(i,j,l)= rui(i,j,l) + source*weight(iw,jw,lw) - elseif (direction.eq.2) then - source= rhop/rhof * volp/volE * f_d(2) - rvi(i,j,l)= rvi(i,j,l) + source*weight(iw,jw,lw) - elseif (direction.eq.3) then - source= rhop/rhof * volp/volE * f_d(3) - rwi(i,j,l)= rwi(i,j,l) + source*weight(iw,jw,lw) - endif + do l=lbeg,lend; do i=ibeg,iend; do j=jbeg,jend + iw=i+1-ibeg + jw=j+1-jbeg + lw=l+1-lbeg + if (direction.eq.1) then + source= rhop/rhof * volp/volE * f_d(1) + rui(i,j,l)= rui(i,j,l) + source*weight(iw,jw,lw) + elseif (direction.eq.2) then + source= rhop/rhof * volp/volE * f_d(2) + rvi(i,j,l)= rvi(i,j,l) + source*weight(iw,jw,lw) + elseif (direction.eq.3) then + source= rhop/rhof * volp/volE * f_d(3) + rwi(i,j,l)= rwi(i,j,l) + source*weight(iw,jw,lw) + endif + enddo; enddo; enddo -5 enddo 2 enddo - 1 enddo ! hg not sure if I want this @@ -321,14 +307,14 @@ ! rwi_o(i,j,l)=rwi(io) !6 enddo - do 6 i=imin,imax; do 6 j=jmin,jmax; do 6 l=lmin,lmax + do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax fsx(i,j,l)= -rui(i,j,l) fsy(i,j,l)= -rvi(i,j,l) fsz(i,j,l)= -rwi(i,j,l) ! fsx(i,j,l)= 0._pr ! fsy(i,j,l)= 0._pr ! fsz(i,j,l)= 0._pr -6 enddo + enddo; enddo; enddo ! remove particles of another processor, they were only needed for source term m=0 @@ -362,9 +348,7 @@ numcol=0 - do 3 i=1,ii - do 3 j=1,jj - do 3 l=1,ll + do i=1,ii; do j=1,jj; do l=1,ll ! condition I only particles in same cell: do 1 m1=2,npic(i,j,l) do 2 m2=1,(m1-1) @@ -480,7 +464,7 @@ 2 enddo 1 enddo -3 enddo + enddo; enddo; enddo numcol=syncSumI(numcol) if (myid.eq.0) write(*,'(a,i9)',advance='no') & @@ -497,7 +481,7 @@ real(kind=pr),allocatable,dimension(:) :: randomx,randomy,randomz integer :: n,i,j,l,ip,jp,lp,npseed,npseedTot, & syncSumI - character*70 :: filename + character(70) :: filename !> number of particles to be seeded npp=np diff --git a/src/post.f90 b/src/post.f90 index 98b11eb..88746c3 100644 --- a/src/post.f90 +++ b/src/post.f90 @@ -36,7 +36,7 @@ syncAv,syncSum,C01,A01 integer :: i,j,l,m,N01,syncSumI logical :: file_ex - character*70 :: filename2 + character(70) :: filename2 diff --git a/src/pre.f90 b/src/pre.f90 index 02a2eb6..ebf46fd 100644 --- a/src/pre.f90 +++ b/src/pre.f90 @@ -31,8 +31,8 @@ use var use mpi integer :: dimitotTemp - character*8 comp_date - character*10 comp_time + character(8) :: comp_date + character(10) :: comp_time ! read input file @@ -343,16 +343,14 @@ zmin=zf(gc); zmax=zf(ll-gc) - do 7 i=2,ii - do 7 j=2,jj - do 7 l=2,ll + do i=2,ii; do j=2,jj; do l=2,ll vol(i,j,l)=(xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1)) -7 enddo + enddo; enddo; enddo volTot= sum(vol(imin:imax,jmin:jmax,lmin:lmax)) ! compute derivatives of grid spacing for mapping - do 8 i=2,ii-1 - if (i.ge.4.and.i.le.ii-3) then + do 8 i=imin,imax + if ((celltype(i+3,jmin,lmin).ne.wall).and.(celltype(i-2,jmin,lmin).ne.wall)) then if (xc(i).gt.0._pr) then call weno(fluxplu,xc(i-2),xc(i-1),xc(i),xc(i+1),xc(i+2)) call weno(fluxmin,xc(i-3),xc(i-2),xc(i-1),xc(i),xc(i+1)) @@ -374,8 +372,8 @@ dxfdi(i)=xc(i+1)-xc(i) endif 8 enddo - do 9 j=2,jj-1 - if (j.ge.4.and.j.le.jj-3) then + do 9 j=jmin,jmax + if ((celltype(imin,j+3,lmin).ne.wall).and.(celltype(imin,j-2,lmin).ne.wall)) then if (yc(j).gt.0._pr) then call weno(fluxplu,yc(j-2),yc(j-1),yc(j),yc(j+1),yc(j+2)) call weno(fluxmin,yc(j-3),yc(j-2),yc(j-1),yc(j),yc(j+1)) @@ -397,8 +395,8 @@ dyfdj(j)=yc(j+1)-yc(j) endif 9 enddo - do 10 l=2,ll-1 - if (l.ge.4.and.l.le.ll-3) then + do 10 l=lmin,lmax + if ((celltype(imin,jmin,l+3).ne.wall).and.(celltype(imin,jmin,l-2).ne.wall)) then if (zc(l).gt.0._pr) then call weno(fluxplu,zc(l-2),zc(l-1),zc(l),zc(l+1),zc(l+2)) call weno(fluxmin,zc(l-3),zc(l-2),zc(l-1),zc(l),zc(l+1)) @@ -443,9 +441,7 @@ ! utau=0.05_pr m=0 - do 1 i=imin,imax - do 1 j=jmin,jmax - do 1 l=lmin,lmax + do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax m=m+1 ! law of the wall/log law: @@ -464,17 +460,17 @@ ! w(i,j,l)=wplus*utau ! linear: - if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l))) - 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)) - u(i,j,l)=1._pr*ubulk*dist/delta + (random(m)-.5_pr)*alpha*ubulk -! u(i,j,l)=3._pr*ubulk*dist/delta -1 enddo + if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l))) + 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)) + u(i,j,l)=1._pr*ubulk*dist/delta + (random(m)-.5_pr)*alpha*ubulk +! u(i,j,l)=3._pr*ubulk*dist/delta + enddo; enddo; enddo v=0._pr w=0._pr u01=u; v01=v; w01=w - call bcUVW + call bcUVWP call sync(u); call sync(v); call sync(w) return diff --git a/src/pressure.f90 b/src/pressure.f90 index 601e140..ac6694d 100644 --- a/src/pressure.f90 +++ b/src/pressure.f90 @@ -21,52 +21,49 @@ u1=u; v1=v; w1=w - do 1 i=imin,imax+1 - do 1 j=jmin,jmax+1 - do 1 l=lmin,lmax+1 - - if (celltype(i,j,l).eq.wall) cycle - - h2sum=0._pr - if (celltype(i+1,j,l).ne.wall) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2 - if (celltype(i-1,j,l).ne.wall) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2 - if (celltype(i,j+1,l).ne.wall) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2 - if (celltype(i,j-1,l).ne.wall) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2 - if (celltype(i,j,l+1).ne.wall) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2 - if (celltype(i,j,l-1).ne.wall) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2 - - - call mass2(ra,i,j,l) - de = (defect_c(i,j,l)+ra)/h2sum*urfp - dex = de/(xf(i)-xf(i-1)) - dey = de/(yf(j)-yf(j-1)) - dez = de/(zf(l)-zf(l-1)) - - if (celltype(i,j,l).eq.active.and.celltype(i+1,j,l).ne.wall) u1(i,j,l) = u1(i,j,l)+dex - if (celltype(i-1,j,l).eq.active.and.celltype(i,j,l).ne.wall) u1(i-1,j,l) = u1(i-1,j,l)-dex - if (celltype(i,j,l).eq.active.and.celltype(i,j+1,l).ne.wall) v1(i,j,l) = v1(i,j,l)+dey - if (celltype(i,j-1,l).eq.active.and.celltype(i,j,l).ne.wall) v1(i,j-1,l) = v1(i,j-1,l)-dey - if (celltype(i,j,l).eq.active.and.celltype(i,j,l+1).ne.wall) w1(i,j,l) = w1(i,j,l)+dez - if (celltype(i,j,l-1).eq.active.and.celltype(i,j,l).ne.wall) w1(i,j,l-1) = w1(i,j,l-1)-dez - - dpa = de*(1._pr/dt+2._pr*nuf/((xf(i)-xf(i-1))**2+(yf(j)-yf(j-1))**2+(zf(l)-zf(l-1))**2)) - - p(i,j,l)= p(i,j,l)+dpa - - da=dpa - if(dpm.lt.da) dpm=da - erru=erru+dex*dex - errv=errv+dey*dey - errw=errw+dez*dez - err=err+da*da - - - 1 enddo + do i=imin,imax+1; do j=jmin,jmax+1; do l=lmin,lmax+1 + + if (celltype(i,j,l).eq.wall) cycle + + h2sum=0._pr + if (celltype(i+1,j,l).ne.wall) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2 + if (celltype(i-1,j,l).ne.wall) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2 + if (celltype(i,j+1,l).ne.wall) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2 + if (celltype(i,j-1,l).ne.wall) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2 + if (celltype(i,j,l+1).ne.wall) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2 + if (celltype(i,j,l-1).ne.wall) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2 + + + call mass2(ra,i,j,l) + de = (defect_c(i,j,l)+ra)/h2sum*urfp + dex = de/(xf(i)-xf(i-1)) + dey = de/(yf(j)-yf(j-1)) + dez = de/(zf(l)-zf(l-1)) + + if (celltype(i,j,l).eq.active.and.celltype(i+1,j,l).ne.wall) u1(i,j,l) = u1(i,j,l)+dex + if (celltype(i-1,j,l).eq.active.and.celltype(i,j,l).ne.wall) u1(i-1,j,l) = u1(i-1,j,l)-dex + if (celltype(i,j,l).eq.active.and.celltype(i,j+1,l).ne.wall) v1(i,j,l) = v1(i,j,l)+dey + if (celltype(i,j-1,l).eq.active.and.celltype(i,j,l).ne.wall) v1(i,j-1,l) = v1(i,j-1,l)-dey + if (celltype(i,j,l).eq.active.and.celltype(i,j,l+1).ne.wall) w1(i,j,l) = w1(i,j,l)+dez + if (celltype(i,j,l-1).eq.active.and.celltype(i,j,l).ne.wall) w1(i,j,l-1) = w1(i,j,l-1)-dez + + dpa = de*(1._pr/dt+2._pr*nuf/((xf(i)-xf(i-1))**2+(yf(j)-yf(j-1))**2+(zf(l)-zf(l-1))**2)) + + p(i,j,l)= p(i,j,l)+dpa + + da=dpa + if(dpm.lt.da) dpm=da + erru=erru+dex*dex + errv=errv+dey*dey + errw=errw+dez*dez + err=err+da*da + + enddo; enddo; enddo u=u1; v=v1; w=w1 call sync(u); call sync(v); call sync(w); call sync(p) - call bcUVW + call bcUVWP err=(syncSum(err)/dimgptot)**(0.5_pr) ! write(*,'(x,a,4(es9.2e2))') 'res. inner it. pressure corr. =',err diff --git a/src/restart.f90 b/src/restart.f90 index 07e41eb..ad3142f 100644 --- a/src/restart.f90 +++ b/src/restart.f90 @@ -22,7 +22,7 @@ subroutine saveField use var integer :: i,j,l,n - character*70 filename + character(70) :: filename write(filename,'(a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,i6.6)') & 'restart/fluidField_p',myid,'_',dimi,'_',dimj,'_',diml,'_',nt @@ -47,7 +47,7 @@ subroutine readField use var integer :: i,j,l,n - character*70 :: filename + character(70) :: filename if (myid.eq.0) write(*,'(a)') 'read old data' diff --git a/src/var.f90 b/src/var.f90 index d660a18..4120fc6 100644 --- a/src/var.f90 +++ b/src/var.f90 @@ -29,7 +29,7 @@ nyw=0.28_pr, & !< duct Poisson ratio Qaccfactor= 10._pr !< artificially accelerate the charging rate - character*70 :: version='pafiX v1.1.0 (Copyright 2019 by H. Grosshans)' + character(70) :: version='pafiX v1.1.0 (Copyright 2019 by H. Grosshans)' ! input real(kind=pr) :: & @@ -74,7 +74,7 @@ ntend !< last time step of computation ! grid - character*1 :: & + character(1) :: & bcx,bcy,bcz, & !< type of boundary condition [(w)all/(p)eriodic] gridx,gridy,gridz !< grid [(u)niform/(s)tretch] real(kind=pr), allocatable, dimension(:,:,:) :: & diff --git a/src/writedat_fluid.f90 b/src/writedat_fluid.f90 index ad191a4..7798953 100644 --- a/src/writedat_fluid.f90 +++ b/src/writedat_fluid.f90 @@ -2,7 +2,7 @@ !> @author Holger Grosshans subroutine writedat_ufluid_xz use var - character*70 filename,rowfmt + character(70) :: filename,rowfmt integer :: i,j,l 100 format(es11.4e2) @@ -33,7 +33,7 @@ !> @author Holger Grosshans subroutine writedat_vfluid_xz use var - character*70 filename,rowfmt + character(70) :: filename,rowfmt integer :: i,j,l 100 format(es11.4e2) @@ -64,7 +64,7 @@ !> @author Holger Grosshans subroutine writedat_wfluid_xz use var - character*70 filename,rowfmt + character(70) :: filename,rowfmt integer :: i,j,l 100 format(es11.4e2) diff --git a/src/writedat_particles.f90 b/src/writedat_particles.f90 index 091bacd..187b225 100644 --- a/src/writedat_particles.f90 +++ b/src/writedat_particles.f90 @@ -5,7 +5,7 @@ use var real(kind=pr) :: partmass integer n - character*80 filename + character(80) :: filename 100 format(10(es13.4e2,x),2(i6,x)) diff --git a/src/writevtk_fluid_xyz.f90 b/src/writevtk_fluid_xyz.f90 index ff8645b..00ee61a 100644 --- a/src/writevtk_fluid_xyz.f90 +++ b/src/writevtk_fluid_xyz.f90 @@ -3,7 +3,7 @@ use var real(kind=pr) :: pp integer :: i,j,l,m - character*70 :: filename,filename2,rowfmt,vecfmt + character(70) :: filename,filename2,rowfmt,vecfmt logical :: file_ex diff --git a/src/writevtk_fluid_xz.f90 b/src/writevtk_fluid_xz.f90 index fbf8130..9c81d44 100644 --- a/src/writevtk_fluid_xz.f90 +++ b/src/writevtk_fluid_xz.f90 @@ -3,7 +3,7 @@ use var real(kind=pr) :: pp integer :: i,j,l,m - character*70 :: filename,filename2,rowfmt,vecfmt + character(70) :: filename,filename2,rowfmt,vecfmt logical :: file_ex diff --git a/src/writevtk_fluid_yz.f90 b/src/writevtk_fluid_yz.f90 index d6130d7..50052ae 100644 --- a/src/writevtk_fluid_yz.f90 +++ b/src/writevtk_fluid_yz.f90 @@ -4,7 +4,7 @@ real(kind=pr), dimension(jj) :: temp4 real(kind=pr) :: temp,temp1,temp2,temp3,pp integer :: i,j,l,m - character*70 :: filename,filename2,rowfmt + character(70) :: filename,filename2,rowfmt logical :: file_ex diff --git a/src/writevtk_grid.f90 b/src/writevtk_grid.f90 index b822d33..346d0e6 100644 --- a/src/writevtk_grid.f90 +++ b/src/writevtk_grid.f90 @@ -2,7 +2,7 @@ subroutine writevtk_grid use var integer :: i,j,l - character*70 filename,rowfmt + character(70) :: filename,rowfmt 100 format(es11.4e2) write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))' diff --git a/src/writevtk_particles.f90 b/src/writevtk_particles.f90 index cc7076c..b61b9d1 100644 --- a/src/writevtk_particles.f90 +++ b/src/writevtk_particles.f90 @@ -4,7 +4,7 @@ subroutine writevtk_particles use var integer :: i,j,l,n,m - character*80 :: filename,filename2,rowfmt,rowfm2 + character(80) :: filename,filename2,rowfmt,rowfm2 logical :: file_ex -- GitLab