diff --git a/src/makefile b/src/makefile index c2b1f4eeab0d4494ec5b7eec11aa0e075d3326f4..7537dc190ebfc4c52b821b740118b81b4892f7f7 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/momentum.f90 b/src/momentum.f90 index c533b544a5e8eb3ef34c6ef2282987c606c65d5e..272ab8c8485296c16d6ff6f0d4e693026de28c44 100644 --- a/src/momentum.f90 +++ b/src/momentum.f90 @@ -94,9 +94,7 @@ integer :: i,j,l - do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax - -! if (celltype(i,j,l).ne.active) cycle + do i=imin,imax+1; do j=jmin,jmax+1; do l=lmin,lmax+1 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. & @@ -108,11 +106,12 @@ defect_c(i,j,l)=0._pr endif + if (celltype(i,j,l).ne.active) cycle + 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 @@ -123,7 +122,6 @@ (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 @@ -134,7 +132,6 @@ (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 diff --git a/src/pre.f90 b/src/pre.f90 index 7b0052160f8d0e2fb64f284e8cbb0a04ee805d70..fd12f704459837b1c138621aeed77d0000e6e734 100644 --- a/src/pre.f90 +++ b/src/pre.f90 @@ -203,6 +203,7 @@ fsx=0._pr; fsy=0._pr; fsz=0._pr defect_c=0._pr; defect_u=0._pr; defect_v=0._pr; defect_w=0._pr celltype=active + dxcdi=0._pr; dycdj=0._pr; dzcdl=0._pr; dxfdi=0._pr; dyfdj=0._pr; dzfdl=0._pr t=0._pr ! electrostatics @@ -350,79 +351,64 @@ ! compute derivatives of grid spacing for mapping do 8 i=imin,imax - if ((celltype(i+2,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)) - else - call weno(fluxplu,xc(i+3),xc(i+2),xc(i+1),xc(i),xc(i-1)) - call weno(fluxmin,xc(i+2),xc(i+1),xc(i),xc(i-1),xc(i-2)) - endif - dxcdi(i)=fluxplu-fluxmin - if (xf(i).gt.0._pr) then - call weno(fluxplu,xf(i-2),xf(i-1),xf(i),xf(i+1),xf(i+2)) - call weno(fluxmin,xf(i-3),xf(i-2),xf(i-1),xf(i),xf(i+1)) - else - call weno(fluxplu,xf(i+3),xf(i+2),xf(i+1),xf(i),xf(i-1)) - call weno(fluxmin,xf(i+2),xf(i+1),xf(i),xf(i-1),xf(i-2)) + if ((celltype(i+2,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)) + else + call weno(fluxplu,xc(i+3),xc(i+2),xc(i+1),xc(i),xc(i-1)) + call weno(fluxmin,xc(i+2),xc(i+1),xc(i),xc(i-1),xc(i-2)) + endif + dxcdi(i)=fluxplu-fluxmin + if (xf(i).gt.0._pr) then + call weno(fluxplu,xf(i-2),xf(i-1),xf(i),xf(i+1),xf(i+2)) + call weno(fluxmin,xf(i-3),xf(i-2),xf(i-1),xf(i),xf(i+1)) + else + call weno(fluxplu,xf(i+3),xf(i+2),xf(i+1),xf(i),xf(i-1)) + call weno(fluxmin,xf(i+2),xf(i+1),xf(i),xf(i-1),xf(i-2)) + endif + dxfdi(i)=fluxplu-fluxmin endif - dxfdi(i)=fluxplu-fluxmin - else - dxcdi(i)=xf(i)-xf(i-1) - dxcdi(i)=xf(i)-xf(i-1) - dxfdi(i)=0._pr - dxfdi(i)=0._pr - endif 8 enddo do 9 j=jmin,jmax - if ((celltype(imin,j+2,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)) - else - call weno(fluxplu,yc(j+3),yc(j+2),yc(j+1),yc(j),yc(j-1)) - call weno(fluxmin,yc(j+2),yc(j+1),yc(j),yc(j-1),yc(j-2)) - endif - dycdj(j)=fluxplu-fluxmin - if (yf(j).gt.0._pr) then - call weno(fluxplu,yf(j-2),yf(j-1),yf(j),yf(j+1),yf(j+2)) - call weno(fluxmin,yf(j-3),yf(j-2),yf(j-1),yf(j),yf(j+1)) - else - call weno(fluxplu,yf(j+3),yf(j+2),yf(j+1),yf(j),yf(j-1)) - call weno(fluxmin,yf(j+2),yf(j+1),yf(j),yf(j-1),yf(j-2)) + if ((celltype(imin,j+2,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)) + else + call weno(fluxplu,yc(j+3),yc(j+2),yc(j+1),yc(j),yc(j-1)) + call weno(fluxmin,yc(j+2),yc(j+1),yc(j),yc(j-1),yc(j-2)) + endif + dycdj(j)=fluxplu-fluxmin + if (yf(j).gt.0._pr) then + call weno(fluxplu,yf(j-2),yf(j-1),yf(j),yf(j+1),yf(j+2)) + call weno(fluxmin,yf(j-3),yf(j-2),yf(j-1),yf(j),yf(j+1)) + else + call weno(fluxplu,yf(j+3),yf(j+2),yf(j+1),yf(j),yf(j-1)) + call weno(fluxmin,yf(j+2),yf(j+1),yf(j),yf(j-1),yf(j-2)) + endif + dyfdj(j)=fluxplu-fluxmin endif - dyfdj(j)=fluxplu-fluxmin - else - dycdj(j)=yf(j)-yf(j-1) - dycdj(j)=yf(j)-yf(j-1) - dyfdj(j)=0._pr - dyfdj(j)=0._pr - endif 9 enddo do 10 l=lmin,lmax - if ((celltype(imin,jmin,l+2).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)) - else - call weno(fluxplu,zc(l+3),zc(l+2),zc(l+1),zc(l),zc(l-1)) - call weno(fluxmin,zc(l+2),zc(l+1),zc(l),zc(l-1),zc(l-2)) - endif - dzcdl(l)=fluxplu-fluxmin - if (zf(l).gt.0._pr) then - call weno(fluxplu,zf(l-2),zf(l-1),zf(l),zf(l+1),zf(l+2)) - call weno(fluxmin,zf(l-3),zf(l-2),zf(l-1),zf(l),zf(l+1)) - else - call weno(fluxplu,zf(l+3),zf(l+2),zf(l+1),zf(l),zf(l-1)) - call weno(fluxmin,zf(l+2),zf(l+1),zf(l),zf(l-1),zf(l-2)) + if ((celltype(imin,jmin,l+2).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)) + else + call weno(fluxplu,zc(l+3),zc(l+2),zc(l+1),zc(l),zc(l-1)) + call weno(fluxmin,zc(l+2),zc(l+1),zc(l),zc(l-1),zc(l-2)) + endif + dzcdl(l)=fluxplu-fluxmin + if (zf(l).gt.0._pr) then + call weno(fluxplu,zf(l-2),zf(l-1),zf(l),zf(l+1),zf(l+2)) + call weno(fluxmin,zf(l-3),zf(l-2),zf(l-1),zf(l),zf(l+1)) + else + call weno(fluxplu,zf(l+3),zf(l+2),zf(l+1),zf(l),zf(l-1)) + call weno(fluxmin,zf(l+2),zf(l+1),zf(l),zf(l-1),zf(l-2)) + endif + dzfdl(l)=fluxplu-fluxmin endif - dzfdl(l)=fluxplu-fluxmin - else - dzcdl(l)=zf(l)-zf(l-1) - dzfdl(l)=zc(l+1)-zc(l) - dzcdl(l)=0._pr - dzfdl(l)=0._pr - endif 10 enddo diff --git a/src/pressure.f90 b/src/pressure.f90 index ac6694da159cfc2aafee1032da72de9e503bf130..37f891568fae22335e341a147682f7e4d57cb908 100644 --- a/src/pressure.f90 +++ b/src/pressure.f90 @@ -1,16 +1,13 @@ !> @author Holger Grosshans !> @brief one sweep to relax mass conservation using distributive !> Jacobi method based on distributive Gauss-Seidel by Brandt (1972) -!> @param direction 1: sweep in increasing index order -!> -1: sweep in decreasing index order !> @param err L2 error norm subroutine pressure(err) use var real(kind=pr),dimension(ii,jj,ll) :: u1,v1,w1 - real(kind=pr) :: err,h2sum,de,dex,dey,dez,ra,dpa,zz,da,dpm, & - erru,errv,errw,syncSum, & - Aw,Ap,Cp,rCp,dpa1,ua,va,wa - integer :: i,j,l,n,nzz + real(kind=pr) :: err,h2sum,de,dex,dey,dez,ra,dpa,dpm, & + erru,errv,errw,syncSum + integer :: i,j,l err=0._pr @@ -32,31 +29,27 @@ 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)) + 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)) + 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 + if (celltype(i,j,l).eq.active) p(i,j,l)= p(i,j,l)+dpa + 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 + if(dpm.lt.dpa) dpm=dpa erru=erru+dex*dex errv=errv+dey*dey errw=errw+dez*dez - err=err+da*da + err=err+dpa*dpa enddo; enddo; enddo