diff --git a/input/input_example.dat b/input/input_example.dat index f8e75a60576e04f220471500628378b1f8e81a17..0aa6a278b383b55b784d8c0bd665ca4a06e1d491 100644 --- a/input/input_example.dat +++ b/input/input_example.dat @@ -14,7 +14,7 @@ interval of writing result files, restart files (-,-) = 10,10000 CFL number (-) = 0.4 -initial bulk flow (m/s), pressure gradient (N/m**3) = +in/initial flow (m/s), pressure gradient (N/m**3) = 3.65,-2.9 fluid density (kg/m**3), fluid kinematic viscosity (m**2/s) = 1.2,1.46E-5 diff --git a/src/bc.f90 b/src/bc.f90 index 091ea2e31536efbb79d476716fba896eab6b1673..d3235856fcc85a844f6019150e5de54da25dc920 100644 --- a/src/bc.f90 +++ b/src/bc.f90 @@ -40,27 +40,13 @@ elseif (bcx.eq.'p') then ! done through the sync routines elseif (bcx.eq.'i') then - if (myid.eq.0) then - call random_number(random) - alpha=.5_pr - mm=0 - 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 + ! see inflow routine if (myid.eq.nrprocs-1) then do 8 m=1,gc - u(imax+m,:,:) =u(imax,:,:) - v(imax+m,:,:) =v(imax,:,:) - w(imax+m,:,:) =w(imax,:,:) + u(imax+m-1,:,:)=u(imax-1,:,:) + v(imax+m,:,:)= v(imax,:,:) + w(imax+m,:,:)= w(imax,:,:) + p(imax+m,:,:)= p(imax,:,:) 8 enddo endif endif @@ -113,6 +99,35 @@ return end +!#################################################################### +!> @author Holger Grosshans +!> @brief inflow boundary condition + subroutine inflow + use var + real(kind=pr),dimension(dimj*diml) :: random + real(kind=pr) :: & + alpha,dist + integer :: i,j,l,m,mm + + + call random_number(random) + alpha=.1_pr + mm=0 + 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)=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 + + return + end + !#################################################################### !> @author Holger Grosshans diff --git a/src/electrostatics.f90 b/src/electrostatics.f90 index 47a3fb31a8bbb5311526d783bb57d7bbdfe63a63..39928c36329e68a6720dea8c50709ef38d2c542c 100644 --- a/src/electrostatics.f90 +++ b/src/electrostatics.f90 @@ -291,8 +291,11 @@ ! pipe: ! ut2=up(n)**2 ! un2=vp(n)**2+wp(n)**2 - - if (direction.eq.2) then +! + if (direction.eq.1) then + ut2=vp(n)**2+wp(n)**2 + un2=up(n)**2 + elseif (direction.eq.2) then ut2=up(n)**2+wp(n)**2 un2=vp(n)**2 elseif (direction.eq.3) then @@ -301,9 +304,11 @@ endif ! John (1979) - alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*(un2+ut2)* & - ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr) - AoAtot=alpha1/4._pr +! 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)) q_el(n)=min(dqp+q_el(n),qpmax) @@ -323,26 +328,27 @@ !> @brief compute particle-particle charging subroutine chargeParticleParticle(n1,n2) use var - real(kind=pr) uax2,urad2,alpha1,dqp + real(kind=pr) urel2,alpha1,dqp integer :: n1,n2 character(70) :: filename - uax2=(up(n1)-up(n2))**2 - urad2=(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2 - - if(urad2.eq.0.) then - dqp=0.5*(q_el(n2)-q_el(n1)) - goto 10 - endif +! urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2 +! +! if (urel2.eq.0._pr) then +! dqp=0.5_pr*(q_el(n2)-q_el(n1)) +! else +!! Soo (1971) +! alpha1=radp(n1)*radp(n2)*(0.625_pr*pi*rhop*(1._pr+restRatio)*urel2 & +! *(radp(n1)+radp(n2))**(0.5_pr)/(radp(n1)**3+radp(n2)**3) & +! *(1._pr-nyp**2)/Ep)**(0.4_pr) +! dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1)) +! endif -! Soo (1971) - alpha1=radp(n1)*radp(n2)*(0.625*pi*rhop*(1+restratio)*(uax2+urad2) & - *(radp(n1)+radp(n2))**(0.5)/(radp(n1)**3+radp(n2)**3) & - *(1-nyp**(2.0))/Ep)**(0.4) + AoAtot=1._pr - dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1)) + dqp=Qaccfactor*AoAtot*(q_el(n2)-q_el(n1)) -10 q_el(n1)=q_el(n1)+dqp + q_el(n1)=q_el(n1)+dqp q_el(n2)=q_el(n2)-dqp diff --git a/src/fluid.f90 b/src/fluid.f90 index 669ea26124f876b16e77832eee005b44674d84ef..b96194da4ec6bb371ce7436e417f2dcec4f4f04b 100644 --- a/src/fluid.f90 +++ b/src/fluid.f90 @@ -9,6 +9,7 @@ call ddt call deferredCorrection + if ((bcx.eq.'i').and.(myid.eq.0)) call inflow do 1 it=1,itmax call momentum(err1) diff --git a/src/makefile b/src/makefile index 7537dc190ebfc4c52b821b740118b81b4892f7f7..6ed0367b799943482043167ac24eec15e7df1821 100644 --- a/src/makefile +++ b/src/makefile @@ -3,21 +3,22 @@ OBJ = var.o main.o pre.o bc.o restart.o fluid.o \ momentum.o post.o schemes.o \ pressure.o writevtk_fluid_xz.o \ - writevtk_fluid_yz.o writedat_fluid.o \ parallel.o writevtk_grid.o particles.o \ mom1.o mom5.o particlesTransport.o \ mass.o electrostatics.o \ writevtk_particles.o writedat_particles.o \ - writevtk_fluid_xyz.o + writevtk_fluid_xy.o writevtk_fluid_xyz.o \ + writevtk_fluid_yz.o writedat_fluid_xz.o \ + writedat_fluid_xy.o 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) @@ -43,14 +44,18 @@ pressure.o: pressure.f90 $(INC) $(CMP) $(FFLAGS) pressure.f90 mass.o: mass.f90 $(INC) $(CMP) $(FFLAGS) mass.f90 +writevtk_fluid_xy.o: writevtk_fluid_xy.f90 $(INC) + $(CMP) $(FFLAGS) writevtk_fluid_xy.f90 writevtk_fluid_xz.o: writevtk_fluid_xz.f90 $(INC) $(CMP) $(FFLAGS) writevtk_fluid_xz.f90 writevtk_fluid_yz.o: writevtk_fluid_yz.f90 $(INC) $(CMP) $(FFLAGS) writevtk_fluid_yz.f90 writevtk_fluid_xyz.o: writevtk_fluid_xyz.f90 $(INC) $(CMP) $(FFLAGS) writevtk_fluid_xyz.f90 -writedat_fluid.o: writedat_fluid.f90 $(INC) - $(CMP) $(FFLAGS) writedat_fluid.f90 +writedat_fluid_xy.o: writedat_fluid_xy.f90 $(INC) + $(CMP) $(FFLAGS) writedat_fluid_xy.f90 +writedat_fluid_xz.o: writedat_fluid_xz.f90 $(INC) + $(CMP) $(FFLAGS) writedat_fluid_xz.f90 writevtk_particles.o: writevtk_particles.f90 $(INC) $(CMP) $(FFLAGS) writevtk_particles.f90 writedat_particles.o: writedat_particles.f90 $(INC) diff --git a/src/momentum.f90 b/src/momentum.f90 index 272ab8c8485296c16d6ff6f0d4e693026de28c44..b1e95c167ea07ca90f369ab586d0d4fbacb54c6d 100644 --- a/src/momentum.f90 +++ b/src/momentum.f90 @@ -13,7 +13,7 @@ integer :: i,j,l err=0._pr - u1=0._pr; v1=0._pr; w1=0._pr + u1=u; v1=v; w1=w dum=0._pr; dvm=0._pr; dwm=0._pr @@ -46,7 +46,6 @@ dwa=abs(dw) if (dwm.lt.dwa) dwm=dwa endif - err=err+dua*dua+dva*dva+dwa*dwa @@ -106,7 +105,7 @@ defect_c(i,j,l)=0._pr endif - if (celltype(i,j,l).ne.active) cycle + if (celltype(i,j,l).ne.active) cycle ! xmax+1 cells are only needed for defect_c 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. & diff --git a/src/parallel.f90 b/src/parallel.f90 index 512fafee15951a0d0989845cb0e6e1c6fb79963c..32bf8753c3c3d09cf1cd84fbde663f87c690f197 100644 --- a/src/parallel.f90 +++ b/src/parallel.f90 @@ -134,6 +134,8 @@ call mpi_wait(rrleft,mpistatus,mpierr) call mpi_wait(rsright,mpistatus,mpierr) call mpi_wait(rrright,mpistatus,mpierr) + if (next.eq.mpi_proc_null) n_recvr=0 + if (prev.eq.mpi_proc_null) n_recvl=0 if (n_sendl.gt.0) call mpi_isend& diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90 index 7619f49af018c710b8fadc3e7f24762bd615acaf..6d274282a7d6f78f0a2005161d133e4083f5c193 100644 --- a/src/particlesTransport.f90 +++ b/src/particlesTransport.f90 @@ -130,7 +130,7 @@ ((myid.eq.nrprocs-1).and.(xp(n).gt.xmax-radp(n)))) then up(n)=-up(n)*restRatio xp(n)=xpold+up(n)*dt - if (qpmax.gt.qp0) call chargeParticleWall(n,1) + if (qpmax.ne.qp0) call chargeParticleWall(n,1) endif elseif (bcx.eq.'i') then if (myid.eq.nrprocs-1) then @@ -150,7 +150,7 @@ vp(n)=-vp(n)*restRatio yp(n)=ypold+vp(n)*dt wcollnum(n)=wcollnum(n)+1 - if (qpmax.gt.qp0) call chargeParticleWall(n,2) + if (qpmax.ne.qp0) call chargeParticleWall(n,2) endif endif @@ -162,7 +162,7 @@ wp(n)=-wp(n)*restRatio zp(n)=zpold+wp(n)*dt wcollnum(n)=wcollnum(n)+1 - if (qpmax.gt.qp0) call chargeParticleWall(n,3) + if (qpmax.ne.qp0) call chargeParticleWall(n,3) endif endif @@ -393,7 +393,7 @@ ppcollnum(n1)=ppcollnum(n1)+1 ppcollnum(n2)=ppcollnum(n2)+1 numcol=numcol+1 -! if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2) + if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2) ! fictitious contact point xpc(n1)=xp(n1)+fracdt*up(n1)*dt diff --git a/src/post.f90 b/src/post.f90 index 88746c3e85f8e20b6dbde9e47f86932f15cabd8f..515f085ff35d92529f748bf77858cc3941d8c478 100644 --- a/src/post.f90 +++ b/src/post.f90 @@ -6,11 +6,15 @@ if (mod(nt,int_results).eq.0.or.(nt).eq.1) then call monitor + call writevtk_fluid_xy call writevtk_fluid_xz call writevtk_fluid_yz - call writedat_ufluid_xz - call writedat_vfluid_xz - call writedat_wfluid_xz + call writedat_ufluid_xy + call writedat_vfluid_xy + call writedat_wfluid_xy +! call writedat_ufluid_xz +! call writedat_vfluid_xz +! call writedat_wfluid_xz ! call writedat_coulomb ! call writevtk_fluid_xyz call writevtk_particles diff --git a/src/pre.f90 b/src/pre.f90 index fd12f704459837b1c138621aeed77d0000e6e734..d95bb5e98d930cfb35fa3a23c2065dec2d82b99b 100644 --- a/src/pre.f90 +++ b/src/pre.f90 @@ -133,7 +133,7 @@ write(20,'(a,es11.2e2)') & 'CFL number (-) = ', cfl write(20,'(a,2(es11.2e2))') & - 'initial bulk flow (m/s), pressure gradient (N/m**3) = ', ubulk,dpdx + 'in/initial flow (m/s), pressure gradient (N/m**3) = ', ubulk,dpdx write(20,'(a,2(es11.2e2))') & 'fluid density (kg/m**3), kinematic visc. (m**2/s) = ', & rhof,nuf @@ -455,8 +455,7 @@ 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 + u(i,j,l)=ubulk*dist/delta + (random(m)-.5_pr)*alpha*ubulk enddo; enddo; enddo v=0._pr diff --git a/src/pressure.f90 b/src/pressure.f90 index 37f891568fae22335e341a147682f7e4d57cb908..a873402fadfc5955bd498d46a99753d19ee3f476 100644 --- a/src/pressure.f90 +++ b/src/pressure.f90 @@ -29,7 +29,7 @@ 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 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)) diff --git a/src/var.f90 b/src/var.f90 index 4120fc64672edfc099523a1d6b6fdefe30f73985..ee6dd1c8a81f372d7947f0dc639777dd7765c488 100644 --- a/src/var.f90 +++ b/src/var.f90 @@ -27,7 +27,7 @@ Ep=1.e8_pr, & !< particle Young's modulus (kg/s**2/m) nyp=0.4_pr , & !< particle Poisson ratio nyw=0.28_pr, & !< duct Poisson ratio - Qaccfactor= 10._pr !< artificially accelerate the charging rate + Qaccfactor= 0.1_pr !< artificially accelerate the charging rate character(70) :: version='pafiX v1.1.0 (Copyright 2019 by H. Grosshans)' diff --git a/src/writedat_fluid_xy.f90 b/src/writedat_fluid_xy.f90 new file mode 100644 index 0000000000000000000000000000000000000000..a3c0e89fbc509af679ccb1756eec5beac996be5c --- /dev/null +++ b/src/writedat_fluid_xy.f90 @@ -0,0 +1,93 @@ +!#################################################################### +!> @author Holger Grosshans + subroutine writedat_ufluid_xy + use var + character(70) :: filename,rowfmt + integer :: i,j,l + +100 format(es11.4e2) + write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))' + + l=int(ll/2) + + write(filename,'(a,i3.3,a,i6.6,a)') & + 'results/fluid_u_xy_p',myid,'_',nt,'.dat' + open(10,file=filename) + + write(10,'(a,es14.6e2,a,x,a)') '# u: t=',t,' s --',version + + write(10,'(a11)',advance='no') '# y/x' + write(10,fmt=rowfmt) (xf(i),i=imin,imax+1) + do j=jmin-1,jmax+1 + write(10,fmt=100,advance='no') yc(j) + write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1) + enddo + + close(10) + + return + + end + +!#################################################################### +!> @author Holger Grosshans + subroutine writedat_vfluid_xy + use var + character(70) :: filename,rowfmt + integer :: i,j,l + +100 format(es11.4e2) + write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))' + + l=int(ll/2) + + write(filename,'(a,i3.3,a,i6.6,a)') & + 'results/fluid_v_xy_p',myid,'_',nt,'.dat' + open(10,file=filename) + + write(10,'(a,es14.6e2,a,x,a)') '# v: t=',t,' s --',version + + write(10,'(a11)',advance='no') '# y/x' + write(10,fmt=rowfmt) (xc(i),i=imin,imax+1) + do j=jmin-1,jmax+1 + write(10,fmt=100,advance='no') yf(j) + write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1) + enddo + + close(10) + + return + + end + +!#################################################################### +!> @author Holger Grosshans + subroutine writedat_wfluid_xy + use var + character(70) :: filename,rowfmt + integer :: i,j,l + +100 format(es11.4e2) + write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))' + + l=int(ll/2) + + write(filename,'(a,i3.3,a,i6.6,a)') & + 'results/fluid_w_xy_p',myid,'_',nt,'.dat' + open(10,file=filename) + + write(10,'(a,es14.6e2,a,x,a)') '# w: t=',t,' s --',version + + write(10,'(a11)',advance='no') '# y/x' + write(10,fmt=rowfmt) (xc(i),i=imin,imax+1) + do j=jmin-1,jmax+1 + write(10,fmt=100,advance='no') yc(j) + write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1) + enddo + + close(10) + + return + + end + diff --git a/src/writedat_fluid.f90 b/src/writedat_fluid_xz.f90 similarity index 100% rename from src/writedat_fluid.f90 rename to src/writedat_fluid_xz.f90 diff --git a/src/writevtk_fluid_xy.f90 b/src/writevtk_fluid_xy.f90 new file mode 100644 index 0000000000000000000000000000000000000000..3ef111e191837c54e130d8b3eddff940733cf2f7 --- /dev/null +++ b/src/writevtk_fluid_xy.f90 @@ -0,0 +1,144 @@ +!> @author Holger Grosshans + subroutine writevtk_fluid_xy + use var + integer :: i,j,l,m + character(70) :: filename,filename2,rowfmt,vecfmt + logical :: file_ex + + +100 format(es11.4e2) + write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))' + write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))' + + l=int(ll/2.)+1 + + write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xy_p',myid,'_',nt,'.vtk' + open(10,file=filename) + +! write visit container file + filename2='results/fluid_xy.visit' + if(myid.eq.0) then + inquire(file=filename2,exist=file_ex) + if (file_ex.and.nt.gt.1) then + open(11,file=filename2,access='append') + else + open(11,file=filename2) + write(11,'(a8,x,i3)') '!NBLOCKS',nrprocs + endif + do m=1,nrprocs + write(11,'(a,i3.3,a,i6.6,a)') 'fluid_xy_p',(m-1),'_',nt,'.vtk' + enddo + close(11) + endif + + + write(10,'(a)') '# vtk DataFile Version 3.0' + write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version + write(10,'(a)') 'ASCII' + write(10,'(a)') 'DATASET RECTILINEAR_GRID' + write(10,'(a10,i4,i4,a4)') & + 'DIMENSIONS',(imax-imin+2),(jmax-jmin+3),'1' + + write(10,*) 'X_COORDINATES',(imax-imin+2),'float' + write(10,fmt=rowfmt) (xc(i),i=imin,imax+1) + + write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float' + write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1) + + write(10,*) 'Z_COORDINATES 1 float' + write(10,100) zc(l) + + write(10,*) 'POINT_DATA ',(imax-imin+2)*(jmax-jmin+3) + + write(10,'(a)') 'SCALARS u float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS v float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS w float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1) + enddo + +! write(10,'(a)') 'SCALARS uraw float 1' +! write(10,'(a)') 'LOOKUP_TABLE default ' +! do j=jmin-1,jmax+1 +! write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1) +! enddo +! +! write(10,'(a)') 'SCALARS vraw float 1' +! write(10,'(a)') 'LOOKUP_TABLE default ' +! do j=jmin-1,jmax+1 +! write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1) +! enddo +! +! write(10,'(a)') 'SCALARS wraw float 1' +! write(10,'(a)') 'LOOKUP_TABLE default ' +! do j=jmin-1,jmax+1 +! write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1) +! enddo + + write(10,'(a)') 'SCALARS p float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (p(i,j,l),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS rho_el float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (rho_el(i,j,l),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS phi_el float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (phi_el(i,j,l),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS Ex_el float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (Ex_el(i,j,l),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS Ey_el float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (Ey_el(i,j,l),i=imin,imax+1) + enddo + + write(10,'(a)') 'SCALARS Ez_el float 1' + write(10,'(a)') 'LOOKUP_TABLE default ' + do j=jmin-1,jmax+1 + write(10,fmt=rowfmt) (Ez_el(i,j,l),i=imin,imax+1) + enddo + +! write(10,'(a)') 'SCALARS Emag_el float 1' +! write(10,'(a)') 'LOOKUP_TABLE default ' +! do j=jmin-1,jmax+1 +! write(10,fmt=rowfmt) & +! (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),i=imin,imax+1) +! enddo +! +! +! write(10,'(a)') 'VECTORS uvw float' +! do j=jmin-1,jmax+1 +! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), & +! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), & +! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1) +! enddo + + + close(10) + + return + end diff --git a/src/writevtk_fluid_xyz.f90 b/src/writevtk_fluid_xyz.f90 index 00ee61ac98674cc6ba610a2deaf87b630fd571d6..22a942959b9325705c337bed4589024a3706244a 100644 --- a/src/writevtk_fluid_xyz.f90 +++ b/src/writevtk_fluid_xyz.f90 @@ -20,7 +20,7 @@ filename2='results/fluid_xyz.visit' if(myid.eq.0) then inquire(file=filename2,exist=file_ex) - if (file_ex.and.nt.ne.0) then + if (file_ex.and.nt.gt.1) then open(11,file=filename2,access='append') else open(11,file=filename2) diff --git a/src/writevtk_fluid_xz.f90 b/src/writevtk_fluid_xz.f90 index 9c81d4434d8dbbf5c7004f177be744d6fc8c829e..832508e46c8367f4b3fb8f6beeafa5f079ae6d4f 100644 --- a/src/writevtk_fluid_xz.f90 +++ b/src/writevtk_fluid_xz.f90 @@ -21,7 +21,7 @@ filename2='results/fluid_xz.visit' if(myid.eq.0) then inquire(file=filename2,exist=file_ex) - if (file_ex.and.nt.ne.0) then + if (file_ex.and.nt.gt.1) then open(11,file=filename2,access='append') else open(11,file=filename2) diff --git a/src/writevtk_fluid_yz.f90 b/src/writevtk_fluid_yz.f90 index 50052aefdba6bd3fa0da711a05fb96ded212f730..59f40d580f4a128435ca9cf0142d611211c78885 100644 --- a/src/writevtk_fluid_yz.f90 +++ b/src/writevtk_fluid_yz.f90 @@ -2,7 +2,7 @@ subroutine writevtk_fluid_yz use var real(kind=pr), dimension(jj) :: temp4 - real(kind=pr) :: temp,temp1,temp2,temp3,pp + real(kind=pr) :: temp,temp1,temp2,temp3 integer :: i,j,l,m character(70) :: filename,filename2,rowfmt logical :: file_ex @@ -12,7 +12,6 @@ write(rowfmt,'(a,i4,a)') '(',(jmax-jmin+3),'(1x,es11.4e2))' i=int(ii/2.)+1 - pp=0._pr write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_yz_p',myid,'_',nt,'.vtk' open(10,file=filename) @@ -21,7 +20,7 @@ filename2='results/fluid_yz.visit' if(myid.eq.0) then inquire(file=filename2,exist=file_ex) - if (file_ex.and.nt.ne.1) then + if (file_ex.and.nt.gt.1) then open(11,file=filename2,access='append') else open(11,file=filename2)