!***************************************************************** ! Generic driving program. Mostly to organize things neatly. ! It does the following things ! ! 1) Call input_data ! 2) Call step ! 3) Call output_data ! 4) Return ! !----------------------------------------------------------------- subroutine drive(pu1,pu2,ptt,paa,pbb,pcc, & cu1,cu2,ctt,caa,cbb,ccc, & ftt,faa,fbb,fcc, & ttt,taa,tbb,tcc, & ucc,vcc, & wcc,xcc,ycc,zcc,& ww1,ww2,ww3,ww4,ww5,ww6,ww7) !----------------------------------------------------------------- use domain_phase use steps_stuff include "mpif.h" real*8, dimension(nl1*nl2+128) :: & & pu1,pu2,ptt,paa,pbb,pcc, & & cu1,cu2,ctt,caa,cbb,ccc, & & ftt,faa,fbb,fcc, & & ttt,taa,tbb,tcc, & & ww1,ww2,ww3,ww4,ww5,ww6,ww7 real*8, dimension(nl1*nl2+128) :: & & ucc,vcc,wcc,xcc,ycc,zcc ! ! Call routine to initialise all the arrays. ! call in_data(pu1,pu2,ptt,paa,pbb,pcc, & & cu1,cu2,ctt,caa,cbb,ccc, & & ftt,faa,fbb,fcc, & & ttt,taa,tbb,tcc, & & wcc,ww1,ww2,ww3,ww4) !write(*,*)'[',nloc,'] done aftr in_data' ! ! Step through the calculation ! do while (nit .lt. ntotal) !write(*,*)'[',nloc,'] calling step' call step(pu1,pu2,ptt,paa,pbb,pcc, & & cu1,cu2,ctt,caa,cbb,ccc, & & ftt,faa,fbb,fcc, & & ttt,taa,tbb,tcc, & & ucc,vcc, & & wcc,xcc,ycc,zcc, & & ww1,ww2,ww3,ww4,ww5,ww6,ww7) !write(*,*)'[',nloc,'] after step' call out_data(pu1,pu2,ptt,paa,pbb,pcc, & & cu1,cu2,ctt,caa,cbb,ccc, & & ww1,ww2,ww3,ww4) end do ! ! Close and leave ! return end subroutine drive !************************************************************** subroutine step(pu1,pu2,ptt,paa,pbb,pcc, & & cu1,cu2,ctt,caa,cbb,ccc, & & ftt,faa,fbb,fcc, & & ttt,taa,tbb,tcc, & & ucc,vcc, & & wcc,xcc,ycc,zcc, & & ww1,ww2,ww3,ww4,ww5,ww6,ww7) !--------------------------------------------------------------- ! ! NOTE order of indices in phase space is (x,y) ! !---------------------------------------------------------------- use domain_phase use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff use output_stuff use gband_stuff use gband_number use communicators use particle use particle_number use input_stuff include 'mpif.h' real*8, dimension(l1s:l1d,l2s:l2d) :: & & pu1,pu2,ptt,paa,pbb,pcc,& & cu1,cu2,ctt,caa,cbb,ccc,& & ftt,faa,fbb,fcc,& & ttt,taa,tbb,tcc,& & ww1,ww2,ww3,ww4,ww5,ww6,ww7 real*8,dimension(l1s:l1d, l2s:l2d) ::& & ucc,vcc,wcc,xcc,ycc,zcc real*8 kx,ky,ks real*8 kxx,kyy,kxy,ksq,cksq real*8 dxx,dyy real*8 status,xx,yy real*8 G, epsilon, A11, A12, A22 real*8 h1,h2,hx,hy,fene integer idir,itype,indi,ierr integer one, root,gi,tip integer shit,atp,i,j,k,wcase integer jj,npt,indicator fene = 2.0d0*la one = 1 root = 0 nit = nit+1 ! write(*,*)'starting step' if(nloc .eq. 0) then if(mod(nit, 10) .eq. 0) write(*,*)'in step',nit,fene,la endif if (falpha .eq. 1) then call velocity1_p(ptt,pu1,pu2,cu1,cu2,ww1,ww2,ww3) else call velocity2_p(ptt,pu1,pu2,cu1,cu2,ww1,ww2,ww3) endif idir = -1; itype = 1 call transform(itype,idir,ptt,ctt,ww1) idir = -1; itype = 1 call transform(itype,idir,pu1,cu1,ww1) idir = -1; itype = 1 call transform(itype,idir,pu2,cu2,ww1) idir = -1; itype = 1 call transform(itype,idir,paa,caa,ww1) call transform(itype,idir,pbb,cbb,ww1) call transform(itype,idir,pcc,ccc,ww1) ! ! Compute the new timestep ! ! call shear_init(ucc,vcc,zcc) if(nit.eq.1) then idir = -1; itype = 1 call transform(itype,idir,ptt,ctt,ww1) if(prev_run_no .eq. 0) then call shear_init(ucc,vcc,zcc) !B-Jet's base flow do i=l1s, l1e do j=l2s, l2e if(vcc(i,j) .ne. 0.0d0) then caa(i,j) = (fene**2.-la**2.+fene**4./4./vcc(i,j)**2.)-& & (fene**4./4./vcc(i,j)**2.)* & & (1.0d0+(8.*fene**2.-16.*la**2.)*vcc(i,j)**2./fene**4.0d0)**0.5 endif if(vcc(i,j) .eq. 0.0d0) then caa(i,j) = la**2.0d0 endif ! write(*,*)'caa(i,j) = ', caa(i,j),i,j enddo enddo ccc = ccc*0.0d0+la**2.0d0 do i=l1s, l1e do j=l2s, l2e cbb(i,j) = la*vcc(i,j)*(fene**2.-(caa(i,j)+ccc(i,j)))/fene**2.0d0 enddo enddo ! caa = 2.*vcc**2.0d0+la**2.0d0 ! cbb = la*vcc idir = 1; itype = 1 call transform(itype,idir,caa,paa,ww1) call transform(itype,idir,cbb,pbb,ww1) call transform(itype,idir,ccc,pcc,ww1) !B-Jet's base flow endif endif call shear_init(ucc,vcc,zcc) ww1 = cu1 + ucc call tmstep(ww1,cu2,caa,cbb,ccc,ww5,ww6) ! if(nloc.eq.0) write(*,*)' after tmstep, dt =',dt !write(*,*)'[',nloc,'] calling diffuse' call diffuse() !write(*,*)'[',nloc,'] done diffuse' ! ! Unload as much of the fluxes as possible ! ww7 = ptt ww1 = paa ww2 = pbb ww3 = pcc do j=l2s,l2e do i=l1s,l1e ptt(i,j) = ptt(i,j) & & + cc1*ftt(i,j)*ext_t1(i,j) & & + cc2*ttt(i,j)*ext_t2(i,j) paa(i,j) = paa(i,j) & & + cc1*faa(i,j)*exs_t1(i,j) & & + cc2*taa(i,j)*exs_t2(i,j) pbb(i,j) = pbb(i,j) & & + cc1*fbb(i,j)*exs_t1(i,j) & & + cc2*tbb(i,j)*exs_t2(i,j) pcc(i,j) = pcc(i,j) & & + cc1*fcc(i,j)*exs_t1(i,j) & & + cc2*tcc(i,j)*exs_t2(i,j) end do end do ! !Swap fluxes around ! ttt = ftt taa = faa tbb = fbb tcc = fcc ! !Start computing the fluxes at the present timestep ! mypi = 4.*atan(1.0)*2. !ucc is U, vcc is U', zcc is omega0 ! call shear_init(ucc,vcc,zcc) ! !Nonlinear Flux ! call shear_init(ucc,vcc,zcc) do j=l2s,l2e do i=l1s,l1e ww4(i,j) = -(cu1(i,j)+ucc(i,j))*caa(i,j) ww5(i,j) = -cu2(i,j)*caa(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww4,wcc,ftt) call transform(itype,idir,ww5,xcc,ftt) do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ww4(i, j) = -kx * wcc(i+1,j) - ky * xcc(i+1,j) ww4(i+1,j) = kx * wcc(i, j) + ky * xcc(i, j) enddo enddo do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ww5(i, j) = -kx * pu1(i+1,j) ww5(i+1,j) = kx * pu1(i, j) ww6(i, j) = -ky * pu1(i+1,j) ww6(i+1,j) = ky * pu1(i, j) enddo enddo idir = -1 ; itype = 1 call transform(itype,idir,ww5,wcc,ftt) call transform(itype,idir,ww6,xcc,ftt) do i=l1s, l1e do j=l2s, l2e ww5(i,j)=2.*caa(i,j)*wcc(i,j)+2.*cbb(i,j)*xcc(i,j)+ & & la**3.0d0*(fene**2./(fene**2.-(caa(i,j)+ccc(i,j))))+ & & la*(1.0d0-fene**2./(fene**2.-(caa(i,j)+ccc(i,j))))*caa(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww5,wcc,ftt) do j=l2s, l2d do i=l1s, l1d faa(i,j) = wcc(i,j)+ww4(i,j) enddo enddo call shear_init(ucc,vcc,zcc) do j=l2s,l2e do i=l1s,l1e ww4(i,j) = -(cu1(i,j)+ucc(i,j))*cbb(i,j) ww5(i,j) = -cu2(i,j)*cbb(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww4,wcc,ftt) call transform(itype,idir,ww5,xcc,ftt) do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ww4(i, j) = -kx * wcc(i+1,j) - ky * xcc(i+1,j) ww4(i+1,j) = kx * wcc(i, j) + ky * xcc(i, j) enddo enddo do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ww5(i, j) = -kx * pu2(i+1,j) ww5(i+1,j) = kx * pu2(i, j) ww6(i, j) = -ky * pu1(i+1,j) ww6(i+1,j) = ky * pu1(i, j) enddo enddo idir = -1 ; itype = 1 call transform(itype,idir,ww5,wcc,ftt) call transform(itype,idir,ww6,xcc,ftt) do i=l1s, l1e do j=l2s, l2e ww5(i,j) = caa(i,j)*wcc(i,j)+ ccc(i,j)*xcc(i,j) + & & la*(1.0d0-fene**2.0d0/(fene**2.-(caa(i,j)+ccc(i,j))))*cbb(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww5,wcc,ftt) do j=l2s,l2d do i=l1s,l1d fbb(i,j) = wcc(i,j)+ww4(i,j) enddo enddo call shear_init(ucc,vcc,zcc) do j=l2s,l2e do i=l1s,l1e ww4(i,j) = -(cu1(i,j)+ucc(i,j))*ccc(i,j) ww5(i,j) = -cu2(i,j)*ccc(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww4,wcc,ftt) call transform(itype,idir,ww5,xcc,ftt) do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ww4(i, j) = -kx * wcc(i+1,j) - ky * xcc(i+1,j) ww4(i+1,j) = kx * wcc(i, j) + ky * xcc(i, j) enddo enddo do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ww5(i, j) = -kx * pu2(i+1,j) ww5(i+1,j) = kx * pu2(i, j) ww6(i, j) = -ky * pu2(i+1,j) ww6(i+1,j) = ky * pu2(i, j) enddo enddo idir = -1 ; itype = 1 call transform(itype,idir,ww5,wcc,ftt) call transform(itype,idir,ww6,xcc,ftt) do i=l1s, l1e do j=l2s, l2e ww5(i,j) = 2.*cbb(i,j)*wcc(i,j)+ & & 2.*ccc(i,j)*xcc(i,j)+ & & la**3.0d0*(fene**2.0d0/(fene**2.-(caa(i,j)+ccc(i,j))))+ & & la*(1.0d0-fene**2.d0/(fene**2.0d0-(caa(i,j)+ccc(i,j))))*ccc(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww5,wcc,ftt) do j=l2s,l2d do i=l1s,l1d fcc(i,j) = ww4(i,j) + wcc(i,j) enddo enddo call shear_init(ucc,vcc,zcc) do j=l2s,l2e do i=l1s,l1e ww4(i,j) = -(cu1(i,j)+ucc(i,j))*(ctt(i,j)+zcc(i,j)) ww5(i,j) = -cu2(i,j)*(ctt(i,j)+zcc(i,j)) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww4,wcc,ftt) call transform(itype,idir,ww5,xcc,ftt) call shear_init(ucc,vcc,zcc) ucc = (0.*1.0d0/Re+ME*la)*ucc call transform(itype,idir,ucc,ww5,ftt) ! ww5 = ww5 * 0.0d0 ftt = ftt*0.0d0 do j=l2s,l2e kx = sk2(j) do i=l1s, l1e-1,2 ky = sk1(i) ftt(i, j) = -kx * wcc(i+1,j) - ky * xcc(i+1,j) & & - ky*(-ky*ky*ww5(i+1,j)) ftt(i+1,j) = kx * wcc(i, j) + ky * xcc(i, j) & & + ky*(-ky*ky*ww5(i, j)) enddo enddo do j=l2s,l2e do i=l1s,l1e ww5(i,j) = (fene**2./(fene**2.-caa(i,j)-ccc(i,j)))*(caa(i,j)-la**2.0d0) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww5,ww1,ww4) do j=l2s,l2e do i=l1s,l1e ww5(i,j) = (fene**2./(fene**2.-caa(i,j)-ccc(i,j)))*cbb(i,j) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww5,ww2,ww4) do j=l2s, l2e do i=l1s, l1e ww5(i,j) = (fene**2./(fene**2.-caa(i,j)-ccc(i,j)))*(ccc(i,j)-la**2.0d0) enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww5,ww3,ww4) do j=l2s,l2e kx = sk2(j) do i=l1s, l1e ky = sk1(i) ftt(i,j) = ftt(i,j) + & & 1.0d0*ME*(kx*ky*(ww1(i,j)-ww3(i,j))+ & & (-kx**2.0d0+ky**2.0d0)*ww2(i,j)) enddo enddo ! !FLUX ! do j=l2s,l2e do i=l1s,l1e ptt(i,j) = (ptt(i,j) & & + cc0*ftt(i,j) )*ext_t0(i,j) paa(i,j) = (paa(i,j) & & + cc0*faa(i,j) )*exs_t0(i,j) pbb(i,j) = (pbb(i,j) & & + cc0*fbb(i,j) )*exs_t0(i,j) pcc(i,j) = (pcc(i,j) & & + cc0*fcc(i,j) )*exs_t0(i,j) end do end do do j=l2s,l2e do i=l1s,l1e if (al_mask(i,j)) then ptt(i,j) = 0.0 paa(i,j) = 0.0 pbb(i,j) = 0.0 pcc(i,j) = 0.0 endif end do end do !test first with all non-Newtonian stresses zero ! paa = 0.0d0 ; pbb = 0.0d0 ; pcc = 0.0d0 ! caa = 0.0d0 ; cbb = 0.0d0 ; ccc = 0.0d0 idir=-1; itype=1 call transform(itype,idir,ptt,ctt,ww1) call transform(itype,idir,paa,caa,ww1) call transform(itype,idir,pbb,cbb,ww1) call transform(itype,idir,pcc,ccc,ww1) mypi = 2.*atan(1.0)*4.0 ww1 = 0.0; ww2 = 0.0; ww3 = 0.0; ww4 = 0.0 ! call shear_init(ucc,vcc,zcc) ! cu1 = ucc ; cu2 = vcc ; ccc=zcc ! Increase time and get out 1000 time = time + dt return end subroutine step subroutine redist_particles() use particle use particle_number use grid_stuff use communicators implicit none include 'mpif.h' integer :: i,ip,ierr,iproc,npart,index integer,allocatable :: my_int_buf(:),int_buf(:) real*8,allocatable :: my_real_buf(:),real_buf(:) integer :: int_buf_size,real_buf_size,np_last_copy integer,allocatable :: recvcounts(:),displs(:) !write(*,*) '[',nloc,'] I have np_last = ',np_last allocate(recvcounts(nproc)) allocate(displs(nproc)) ! error check if ( np_last /= mypp(2) + mymp(2)) then write(*,*) 'Error: np_last /= mypp(2) + mymp(2).' call MPI_FINALIZE(ierr) stop end if ! first, figure out process for each particle... do ip = 1,np_last lparticle(ip)%ip = get_np_for_xy(lparticle(ip)%newx,lparticle(ip)%newy) end do ! cycle through all processes and figure out the max local buf sizes we will need... real_buf_size = 0 int_buf_size = 0 do iproc = 0,nproc-1 if (nloc /= iproc) then ! check how many particles we will be sending to iproc... npart = count(lparticle(1:np_last)%ip == iproc) real_buf_size = max(real_buf_size,N_PASSED_REALS*npart) int_buf_size = max(int_buf_size,N_PASSED_INTS*npart) ! write(*,*) '[',nloc,'] sending ',npart,' out of ',np_last,' to proc ',iproc end if end do ! write(*,*) '[',nloc,'] made it 1, int_buf_size = ',int_buf_size,' real_buf_size = ',real_buf_size !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) allocate(my_int_buf(int_buf_size)) allocate(my_real_buf(real_buf_size)) !write(*,*) '[',nloc,'] made it 2' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! send particles to new process, one process at a time... do iproc = 0,nproc-1 if (nloc /= iproc) then int_buf_size = 0 real_buf_size = 0 npart = 0 do ip = 1,np_last if (lparticle(ip)%ip == iproc) then npart = npart + 1 call pack_p_reals(lparticle(ip), & my_real_buf(real_buf_size+1:real_buf_size+N_PASSED_REALS)) real_buf_size = real_buf_size + N_PASSED_REALS call pack_p_ints(lparticle(ip), & my_int_buf(int_buf_size+1:int_buf_size+N_PASSED_INTS)) int_buf_size = int_buf_size + N_PASSED_INTS end if end do else npart = 0 end if !write(*,*) '[',nloc,'] made it 3' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! gather the sizes at iproc call MPI_GATHER(npart,1,MPI_INTEGER, & recvcounts,1,MPI_INTEGER, & iproc,MPI_COMM_WORLD,ierr) !write(*,*) '[',nloc,'] made it 4' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! reals first... if (nloc == iproc) then !write(*,*) '[',nloc,'] recvcounts = ',recvcounts(1:nproc) ! we are the process that will be getting the particles... npart = sum(recvcounts(1:nproc)) recvcounts(1:nproc) = recvcounts(1:nproc)*N_PASSED_REALS displs(1) = 0 do i = 2, nproc displs(i) = displs(i-1) + recvcounts(i-1) end do real_buf_size = displs(nproc) + recvcounts(nproc) allocate(real_buf(real_buf_size)) real_buf_size = 0 !write(*,*) '[',nloc,'] recvcounts = ',recvcounts(1:nproc) !write(*,*) '[',nloc,'] displs = ',displs(1:nproc) end if !write(*,*) '[',nloc,'] made it 5, real_buf_size = ',real_buf_size !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) call MPI_GATHERV(my_real_buf,real_buf_size,MPI_REAL8, & real_buf,recvcounts,displs,MPI_REAL8, & iproc,MPI_COMM_WORLD,ierr) !write(*,*) '[',nloc,'] made it 6' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) if (nloc == iproc) then ! unpack into NEW particles... ! store np_last... np_last_copy = np_last real_buf_size = 0 do ip = 1,npart index = get_new_p() call unpack_p_reals(real_buf(real_buf_size+1:real_buf_size+N_PASSED_REALS),lparticle(index)) real_buf_size = real_buf_size + N_PASSED_REALS end do end if !write(*,*) '[',nloc,'] made it 7' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! ints next... if (nloc == iproc) then recvcounts(1:nproc) = recvcounts(1:nproc)/N_PASSED_REALS*N_PASSED_INTS displs(1) = 0 do i = 2, nproc displs(i) = displs(i-1) + recvcounts(i-1) end do int_buf_size = displs(nproc) + recvcounts(nproc) allocate(int_buf(int_buf_size)) int_buf_size = 0 end if !write(*,*) '[',nloc,'] made it 8, int_buf_size = ',int_buf_size !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) call MPI_GATHERV(my_int_buf,int_buf_size,MPI_INTEGER, & int_buf,recvcounts,displs,MPI_INTEGER, & iproc,MPI_COMM_WORLD,ierr) !write(*,*) '[',nloc,'] made it 9' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) if (nloc == iproc) then ! unpack into last particles... ! error check... if (np_last_copy + npart /= np_last) then write(*,*) 'Error: (np_last_copy + npart /= np_last).' call MPI_FINALIZE(ierr) stop end if int_buf_size = 0 do index = np_last_copy+1,np_last_copy+npart call unpack_p_ints(int_buf(int_buf_size+1:int_buf_size+N_PASSED_INTS),lparticle(index)) int_buf_size = int_buf_size + N_PASSED_INTS if (lparticle(index)%sign == 1) then mypp(2) = mypp(2) + 1 else mymp(2) = mymp(2) + 1 end if lparticle(index)%ip = nloc end do end if call MPI_BARRIER(MPI_COMM_WORLD,ierr) end do !write(*,*) '[',nloc,'] made it 10' !if (nloc == 0) read(*,*) ierr ! clean up the particle vector... index = 0 do ip = 1,np_last if (lparticle(ip)%ip /= nloc) then if (lparticle(ip)%sign == 1) then mypp(2) = mypp(2) - 1 else mymp(2) = mymp(2) - 1 end if index = index + 1 else if (index > 0) then call copy_p(lparticle(ip),lparticle(ip-index)) end if end do np_last = np_last - index ! check... if (np_last /= mypp(2) + mymp(2)) then write(*,*) 'Error: (np_last_copy + npart /= np_last).' call MPI_FINALIZE(ierr) stop end if !write(*,*) '[',nloc,'] I now have np_last = ',np_last !write(*,*) '[',nloc,'] made it 11' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) if (allocated(my_int_buf)) deallocate(my_int_buf) if (allocated(my_real_buf)) deallocate(my_real_buf) if (allocated(int_buf)) deallocate(int_buf) if (allocated(real_buf)) deallocate(real_buf) deallocate(recvcounts) deallocate(displs) !write(*,*) '[',nloc,'] made it 12' !if (nloc == 0) read(*,*) ierr !call MPI_BARRIER(MPI_COMM_WORLD,ierr) end subroutine redist_particles !************************************************************** subroutine tmstep(cu1,cu2,caa,cbb,ccc,ww1,ww2) !--------------------------------------------------------- use domain_confg use time_stuff use grid_stuff use steps_stuff use physics_stuff real*8, dimension(l1s:l1d,l2s:l2d) :: & & cu1,cu2,ww1,ww2 real*8, dimension(l1s:l1d,l2s:l2d) :: caa,cbb,ccc real*8, dimension(1) :: tmp1,tmp2 real*8, dimension(1) :: tmp3,tmp4 real*8, dimension(1) :: tmp5,tmp6 real*8, dimension(1) :: tmp7,tmp8 real*8 epsilon, dt_new integer nn, idir, itype ! ! If safety > 0.0 the timestep is computed every timestep ! If safety < 0.0 timestep is FIXED and dt= -safety ! if (safety .lt. 0) then dt_new = -safety else ww1 = cu1*cu1 + cu2*cu2 tmp1(1) = maxval(abs(ww1(l1s:l1e,l2s:l2e))) nn = 1 call globalmax(tmp1,nn,tmp2,gcomm) tmp3(1) = 0.0d0 ; tmp5(1) = 0.0d0 ; tmp7(1) = 0.0d0 do i=l1s, l1d do j=l2s, l2d tmp3(1) = tmp3(1) + (caa(i,j)**2.0d0)/ng1/ng2 tmp5(1) = tmp5(1) + (cbb(i,j)**2.0d0)/ng1/ng2 tmp7(1) = tmp7(1) + (ccc(i,j)**2.0d0)/ng1/ng2 enddo enddo nn=1 call globalsum(tmp3,nn,tmp4,gcomm) call globalsum(tmp5,nn,tmp6,gcomm) call globalsum(tmp7,nn,tmp8,gcomm) tmp4(1)=tmp4(1)**0.5 tmp6(1)=tmp6(1)**0.5 tmp8(1)=tmp8(1)**0.5 ! ! tmp2(1) contains the global maximum over all the ! processors of ww1 ! dt_new = safety*min(dx,dy) & /max(0.1,sqrt(tmp2(1))) ! dt_new = min(tmp4(1)/tmp2(1),dt_new) ! dt_new = min(tmp6(1)/tmp2(1),dt_new) ! dt_new = min(tmp8(1)/tmp2(1),dt_new) ! dt_new = min(1.0d0/Re/tmp2(1),dt_new) ! if(nloc.eq.0) write(*,*)'dt_new = ',tmp3,tmp5,tmp7,dt_new end if ! ! Compute time integration constants ! if (nit .eq. 1) then ! 1st step is Euler dt2 = 0.0 dt1 = 0.0 dt = dt_new cc0 = dt cc1 = 0.0 cc2 = 0.0 else if (nit .eq. 2) then ! Second step is 2 levels AB dt2 = 0.0 dt1 = dt dt = dt_new cc0 = dt*(0.5*dt/dt1 + 1.0) cc1 = -0.5*dt*(dt/dt1) cc2 = 0.0 else ! All other steps are 3 levels AB dt2 = dt1 + dt dt1 = dt dt = dt_new cc0 = dt + (dt/dt1)*(dt/dt2)*(dt/3.0 + 0.5*(dt1+ dt2)) cc1 = (dt/dt1)*(dt/(dt1-dt2))*(dt/3.0 + 0.5*dt2) cc2 = -(dt/dt2)*(dt/(dt1-dt2))*(dt/3.0 + 0.5*dt1) end if return end subroutine tmstep !**************************************************************** subroutine diffuse() !---------------------------------------------------------------- use domain_phase use diffuse_stuff use physics_stuff use deriv_stuff use time_stuff real*8 ks, ffactor !-------------------------------------- time level twice removed do i=l1s, l1e do j=l2s, l2e ext_t2(i,j) = ext_t1(i,j) * ext_t0(i,j) exs_t2(i,j) = exs_t1(i,j) * exs_t0(i,j) enddo enddo !-------------------------------------- time level once removed ext_t1 = ext_t0 exs_t1 = exs_t0 !-------------------------------------- present time level !HYPER---> 1012/99 do i=l1s, l1e do j=l2s, l2e ks = ak1(i)**2. + ak2(j)**2. ext_t0(i,j)=exp(-dt*ks/Re) exs_t0(i,j)=exp(-dt*la) enddo enddo ! return end subroutine diffuse !******************************************************************* subroutine limit(ptt,ctt,ww1) !******************************************************************* use domain_confg use fft_stuff ! real*8, dimension(l1s:l1d,l2s:l2d) :: ptt,ctt,ww1 integer idir,itype integer i,j,k ! idir = -1; itype = 1 !(cos) call transform(itype,idir,ptt,ctt,ww1) ! do j=l2s,l2e do i=l1s,l1e ! if (ctt(i,j) .gt. 1.0) then ctt(i,j)=1.0 else if(ctt(i,j) .lt. 0.0) then ctt(i,j)=0.0 end if end do end do ! idir = 1; itype = 1 !(cos) call transform(itype,idir,ctt,ptt,ww1) ! return end subroutine limit !********************************************************** subroutine velocity1_p(ptt,pu1,pu2,cu1,cu2,ww1,ww2,ww3) !********************************************************** use domain_phase use wave_stuff use deriv_stuff use physics_stuff use alias_stuff ! real*8, dimension(l1s:l1d,l2s:l2d) :: & & pu1,pu2,ptt, & & cu1,cu2,ww1,ww2,ww3 ! real*8 kx,ky real*8 kxx,kyy,kk,tk1,tk2,hkk do j=l2s,l2e kx=sk2(j) kxx = kx*kx do i=l1s,l1e-1,2 ky=sk1(i) kyy=ky*ky kk=sqrt(kxx+kyy) ! if (kk .ge. tiny) then hkk = 1.0/kk tk1 = ptt(i ,j)*hkk tk2 = ptt(i+1,j)*hkk ! pu1(i ,j) = -ky*tk2 pu1(i+1,j) = ky*tk1 ! pu2(i ,j) = kx*tk2 pu2(i+1,j) = -kx*tk1 ! else pu1(i ,j) = 0.0 pu1(i+1,j) = 0.0 ! pu2(i ,j) = 0.0 pu2(i+1,j) = 0.0 end if end do end do return end subroutine velocity1_p !********************************************************** subroutine velocity2_p(ptt,pu1,pu2,cu1,cu2,ww1,ww2,ww3) !********************************************************** use domain_phase use wave_stuff use deriv_stuff use physics_stuff use alias_stuff ! real*8, dimension(l1s:l1d,l2s:l2d) :: & & pu1,pu2,ptt, & & cu1,cu2,ww1,ww2,ww3 ! real*8 kx,ky real*8 kxx,kyy,kk,tk1,tk2,hkk do j=l2s,l2e kx=sk2(j) kxx = kx*kx do i=l1s,l1e-1,2 ky=sk1(i) kyy=ky*ky kk= (kxx+kyy) ! if (kk .ge. tiny) then hkk = 1.0/kk !tk = STREAMFUNCTION = VORTICITY/(k^2) tk1 = ptt(i ,j)*hkk tk2 = ptt(i+1,j)*hkk !U = D (STREAMFUNCTION)/ DZ pu1(i ,j) = -ky*tk2 pu1(i+1,j) = ky*tk1 !W =-D (STREAMFUNCTION)/ DX pu2(i ,j) = kx*tk2 pu2(i+1,j) = -kx*tk1 ! else pu1(i ,j) = 0.0 pu1(i+1,j) = 0.0 ! pu2(i ,j) = 0.0 pu2(i+1,j) = 0.0 end if end do end do return end subroutine velocity2_p !********************************************************** subroutine velocity_c(ptt,pu1,pu2,cu1,cu2,ww1,ww2,ww3) !********************************************************** use domain_confg use physics_stuff use grid_stuff use alias_stuff ! real*8, dimension(l1s:l1d,l2s:l2d) :: & & pu1,pu2,ptt, & & cu1,cu2,ww1,ww2,ww3 ! integer i,j real*8 x,y, csy , sny do j=l2s,l2e y=float(j-1)*dy sny = sin(twopi*y/Ymx) csy = cos(twopi*y/Ymx) do i=l1s,l1e x=float(i-1)*dx ! cu1(i,j) = cos(twopi*x/Xmx)*sny*Xmx ! cu2(i,j) = -sin(twopi*x/Xmx)*csy*ymx cu1(i,j) = 0.0 cu2(i,j) = sin(twopi*x/Xmx) end do end do return end subroutine velocity_c subroutine ttd(fss, css, ww4, ww1, ww2, ww3) use domain_phase use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff real*8, dimension(l1s:l1d, l2s:l2d) :: fss, css, ww4, ww1, ww2, ww3 real*8 ky, kyy integer idir, itype do j=l2s, l2e do i=l1s, l1e-1, 2 ky = sk1(i) ww1(i,j) = -ky * ww4(i+1,j) ww1(i+1,j) = ky * ww4(i, j) enddo enddo idir = -1; itype = 1 call transform(itype,idir,ww1,ww2,ww3) ! ww2 is theta in configuration space do j=l2s, l2e do i=l1s, l1e ky = sk1(i) ww1(i, j) = -ky*ky * ww4(i, j) enddo enddo idir = -1; itype = 1 call transform(itype,idir,ww1,css,ww3) ! print *,'in ttd' do i=l1s, l1d do j=l2s, l2d ! print *,i,j,css(i,j) ww2(i, j) = -2.*ww2(i,j)**2*css(i,j)/(1.+ww2(i, j)**2)**2. enddo enddo idir = 1; itype = 1 call transform(itype,idir,ww2,ww3,ww1) ww1 = ww3 css = 0.0 return end subroutine ttd subroutine weno_adv(pcc,cu1,cu2, ww1, ww2) use domain_confg use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff real*8, dimension(l1s:l1d, l2s:l2d) :: cu1,cu2,ww1,ww2 real*8, dimension(l1s-3:l1d+3, l2s-3:l2d+3) :: pcc real*8 eps,nu1,nu2,nu3,nu4,nu5,s1,s2,s3,a1,a2,a3,w1,w2,w3 eps = 10.**(-6) do i=l1s, l1d do j=l2s, l2d if(cu1(i,j) .gt. 0) then nu1 = (pcc(i-2,j)-pcc(i-3,j))/dx nu2 = (pcc(i-1,j)-pcc(i-2,j))/dx nu3 = (pcc(i, j)-pcc(i-1,j))/dx nu4 = (pcc(i+1,j)-pcc(i, j))/dx nu5 = (pcc(i+2,j)-pcc(i+1,j))/dx endif if(cu1(i,j) .le. 0) then nu1 = (pcc(i+3,j)-pcc(i+2,j))/dx nu2 = (pcc(i+2,j)-pcc(i+1,j))/dx nu3 = (pcc(i+1,j)-pcc(i, j))/dx nu4 = (pcc(i, j)-pcc(i-1,j))/dx nu5 = (pcc(i-1,j)-pcc(i-2,j))/dx endif s1 = 13./12.*(nu1-2*nu2+ nu3)**2.+ & & 1. / 4.*(nu1-4*nu2+3*nu3)**2. s2 = 13./12.*(nu2-2*nu3+ nu4)**2.+ & & 1. / 4.*(nu2-nu4)**2. s3 = 13./12.*(nu3-2*nu4+ nu5)**2.+ & & 1. / 4.*(3*nu3-4*nu4+nu5)**2. a1 = (1./10.)/(s1+eps)**2. a2 = (6./10.)/(s2+eps)**2. a3 = (3./10.)/(s3+eps)**2. w1 = a1/(a1+a2+a3) w2 = a2/(a1+a2+a3) w3 = a3/(a1+a2+a3) ww1(i,j) = w1*( nu1/3.-7*nu2/6.+11.*nu3/6.)+ & & w2*(-nu2/6.+5.*nu3/6.+nu4/3.)+ & & w3*( nu3/3.+5.*nu4/6.-nu5/6.) enddo enddo do i=l1s, l1d do j=l2s, l2d if(cu2(i,j) .gt. 0) then nu1 = (pcc(i,j-2)-pcc(i,j-3))/dy nu2 = (pcc(i,j-1)-pcc(i,j-2))/dy nu3 = (pcc(i, j)-pcc(i,j-1))/dy nu4 = (pcc(i,j+1)-pcc(i, j))/dy nu5 = (pcc(i,j+2)-pcc(i,j+1))/dy endif if(cu2(i,j) .le. 0) then nu1 = (pcc(i,j+3)-pcc(i,j+2))/dy nu2 = (pcc(i,j+2)-pcc(i,j+1))/dy nu3 = (pcc(i,j+1)-pcc(i, j))/dy nu4 = (pcc(i, j)-pcc(i,j-1))/dy nu5 = (pcc(i,j-1)-pcc(i,j-2))/dy endif s1 = 13./12.*(nu1-2*nu2+ nu3)**2.+ & & 1. / 4.*(nu1-4*nu2+3*nu3)**2. s2 = 13./12.*(nu2-2*nu3+ nu4)**2.+ & & 1. / 4.*(nu2-nu4)**2. s3 = 13./12.*(nu3-2*nu4+ nu5)**2.+ & & 1. / 4.*(3*nu3-4*nu4+nu5)**2. a1 = 1./10./(s1+eps)**2. a2 = 6./10./(s2+eps)**2. a3 = 3./10./(s3+eps)**2. w1 = a1/(a1+a2+a3) w2 = a2/(a1+a2+a3) w3 = a3/(a1+a2+a3) ww2(i,j) = w1*( nu1/3.-7*nu2/6.+11.*nu3/6.)+ & & w2*(-nu2/6.+5.*nu3/6.+nu4/3.)+ & & w3*( nu3/3.+5.*nu4/6.-nu5/6.) enddo enddo return end subroutine weno_adv subroutine mp_level(pcc,ww1,ww2,ww3,ww4) use communicators use domain_confg use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff include 'mpif.h' integer mpcount real*8, dimension(l1s-3:l1d+3, 3) :: ww1,ww2,ww3,ww4 real*8, dimension(l1s-3:l1d+3, l2s-3:l2d+3) :: pcc integer status, ierr integer status_array(MPI_STATUS_SIZE,4), ierr,req(4) do j=l2s, l2d pcc(l1d+1,j) = pcc(l1s+0,j) pcc(l1d+2,j) = pcc(l1s+1,j) pcc(l1d+3,j) = pcc(l1s+2,j) pcc(l1s-1,j) = pcc(l1d+0,j) pcc(l1s-2,j) = pcc(l1d-1,j) pcc(l1s-3,j) = pcc(l1d-2,j) enddo if(nproc.eq.1) then do i=l1s-3, l1d+3 pcc(i,l2d+1) = pcc(i,l2s+0) pcc(i,l2d+2) = pcc(i,l2s+1) pcc(i,l2d+3) = pcc(i,l2s+2) pcc(i,l2s-1) = pcc(i,l2d+0) pcc(i,l2s-2) = pcc(i,l2d-1) pcc(i,l2s-3) = pcc(i,l2d-2) enddo endif mpcount = (l1d-l1s+7)*3 if(mod(nloc,2)==0) then ww1(l1s-3:l1d+3,1:3) = pcc(l1s-3:l1d+3,l2d-2:l2d ) ww3(l1s-3:l1d+3,1:3) = pcc(l1s-3:l1d+3,l2s: l2s+2) else ww4(l1s-3:l1d+3,1:3) = pcc(l1s-3:l1d+3,l2d-2:l2d ) ww2(l1s-3:l1d+3,1:3) = pcc(l1s-3:l1d+3,l2s: l2s+2) endif if(nproc.gt.1) then if(mod(nloc,2)==0) then nloc_nbr = mod(nloc + 1, nproc) call MPI_SEND(ww1,mpcount,MPI_REAL8,nloc_nbr,MY_TAG,& & MPI_COMM_WORLD,ierr) call MPI_RECV(ww2,mpcount,MPI_REAL8,nloc_nbr,& & MY_TAG+1,MPI_COMM_WORLD,status,ierr) nloc_nbr = nloc-1 if(nloc_nbr == -1) nloc_nbr = nproc-1 call MPI_SEND(ww3,mpcount,MPI_REAL8,nloc_nbr,MY_TAG+2,& & MPI_COMM_WORLD,ierr) call MPI_RECV(ww4,mpcount,MPI_REAL8,nloc_nbr,MY_TAG+3, & & MPI_COMM_WORLD,status, ierr) else nloc_nbr = nloc-1 if(nloc_nbr .lt.0) nloc_nbr=nloc_nbr+nproc call MPI_RECV(ww1,mpcount,MPI_REAL8,nloc_nbr,MY_TAG, & & MPI_COMM_WORLD,status,ierr) call MPI_SEND(ww2,mpcount,MPI_REAL8,nloc_nbr,MY_TAG+1,& & MPI_COMM_WORLD,ierr) nloc_nbr = mod(nloc+1, nproc) call MPI_RECV(ww3,mpcount,MPI_REAL8,nloc_nbr,MY_TAG+2, & & MPI_COMM_WORLD,status,ierr) call MPI_SEND(ww4,mpcount,MPI_REAL8,nloc_nbr,MY_TAG+3, & & MPI_COMM_WORLD,ierr) endif ! if (nloc==0) read(*,*) ierr ! call MPI_BARRIER(MPI_COMM_WORLD,ierr) ! print *,'just about done in mp_level' if(mod(nloc,2)==0) then pcc(l1s-3:l1d+3,l2d+1:l2d+3) = ww2(l1s-3:l1d+3,1:3) pcc(l1s-3:l1d+3,l2s-3:l2s-1) = ww4(l1s-3:l1d+3,1:3) else pcc(l1s-3:l1d+3,l2d+1:l2d+3) = ww3(l1s-3:l1d+3,1:3) pcc(l1s-3:l1d+3,l2s-3:l2s-1) = ww1(l1s-3:l1d+3,1:3) endif call MPI_BARRIER(MPI_COMM_WORLD,ierr) endif return end subroutine mp_level subroutine RK3_G(pcc,tcc,cu1,cu2,ww1,ww2,ww3,ww4) use communicators use domain_confg use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff include 'mpif.h' integer mpcount,case real*8, dimension(l1s:l1d, l2s:l2d) :: cu1, cu2 real*8, dimension(l1s:l1d, l2s:l2d) :: ww1, ww2, ww3, ww4 real*8, dimension(l1s-3:l1d+3, l2s-3:l2d+3) :: pcc,tcc real*8 status, ierr !write(*,*)'in RK3_G',dt call MPI_BARRIER(MPI_COMM_WORLD,ierr) call mp_level(pcc,ww1,ww2,ww3,ww4) call MPI_BARRIER(MPI_COMM_WORLD,ierr) call weno_adv(pcc,cu1,cu2,ww1,ww2) do i=l1s, l1d do j=l2s, l2d tcc(i,j) = pcc(i,j) - & & dt*(cu1(i,j)*ww1(i,j)+cu2(i,j)*ww2(i,j)) enddo enddo call MPI_BARRIER(MPI_COMM_WORLD,ierr) call mp_level(tcc,ww1,ww2,ww3,ww4) call MPI_BARRIER(MPI_COMM_WORLD,ierr) call weno_adv(tcc,cu1,cu2,ww1,ww2) do i=l1s, l1d do j=l2s, l2d tcc(i,j) = tcc(i,j) - & & dt*(cu1(i,j)*ww1(i,j)+cu2(i,j)*ww2(i,j)) tcc(i,j) = 3./4.*pcc(i,j)+1./4.*tcc(i,j) enddo enddo call MPI_BARRIER(MPI_COMM_WORLD,ierr) call mp_level(tcc,ww1,ww2,ww3,ww4) call MPI_BARRIER(MPI_COMM_WORLD,ierr) call weno_adv(tcc,cu1,cu2,ww1,ww2) do i=l1s, l1d do j=l2s, l2d tcc(i,j) = tcc(i,j) - & & dt*(cu1(i,j)*ww1(i,j)+cu2(i,j)*ww2(i,j)) tcc(i,j) = 1./3.*pcc(i,j)+2./3.*tcc(i,j) enddo enddo return end subroutine RK3_G subroutine weno_phi(pcc, ww1, ww2, case) use domain_confg use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff real*8, dimension(l1s:l1d, l2s:l2d) :: ww1,ww2 real*8, dimension(l1s-3:l1d+3, l2s-3:l2d+3) :: pcc real*8 eps,ap,am,bp,bm,cp,cm,dp,dm,sp,sm real*8 a,b,c,d,s real*8 nu1,nu2,nu3,nu4,nu5,s1,s2,s3,a1,a2,a3,w1,w2,w3 integer xy,dir,i,j,case eps = 10.**(-6) do i=l1s, l1d do j=l2s, l2d xy = 1; dir = 1 call weno_grad(pcc,i,j,xy,dir,a) xy = 1; dir = -1 call weno_grad(pcc,i,j,xy,dir,b) xy = 2; dir = 1 call weno_grad(pcc,i,j,xy,dir,c) xy = 2; dir = -1 call weno_grad(pcc,i,j,xy,dir,d) ap = max(a,0.0); am = min(a,0.0) bp = max(b,0.0); bm = min(b,0.0) cp = max(c,0.0); cm = min(c,0.0) dp = max(d,0.0); dm = min(d,0.0) if(case.eq.1) then ww1(i,j) = sqrt(max(ap**2.,bm**2.)+max(cp**2.,dm**2.))-1. ww2(i,j) = sqrt(max(am**2.,bp**2.)+max(cm**2.,dp**2.))-1. endif if(case.eq.2) then ww1(i,j) = ((a+b)/2.)**2.+((c+d)/2.)**2. endif if(case.eq.3) then if((a+b)**2.+(c+d)**2.gt.0) then ww1(i,j) = (a+b)/sqrt((a+b)**2.+(c+d)**2.) ww2(i,j) = (c+d)/sqrt((a+b)**2.+(c+d)**2.) else ww1(i,j) = 0 ; ww2(i,j) = 0 endif endif if(case.eq.4) then ww1(i,j) = (a+b)/2. ww2(i,j) = (c+d)/2. endif enddo enddo return end subroutine weno_phi subroutine RK3_phi(pcc,fcc,tcc,ww1,ww2,ww3,ww4) use communicators use domain_confg use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff ! include 'mpif.h' integer mpcount real*8, dimension(l1s:l1d, l2s:l2d) :: & & ww1, ww2, ww3, ww4 real*8, dimension(l1s-3:l1d+3, l2s-3:l2d+3) :: pcc,fcc,tcc real*8 sn,sp,sm,ddx,tol,tmp1(1),tmp2(1),tmp3(1),tol_reinit integer ii,iter,one,root,wcase ddx = min(dx,dy) ! print *,' in RK3_phi ',nloc,ddx,dt ! print *,' in RK3_phi, calling weno_phi ',nloc one = 1 root = 0 tmp1(1) = 10.**8. ; tmp2(1) = 10.**8. ; tmp3(1) = 10.**8. tol_reinit =0.001*min(dx,dy) iter = 20 ii = 0 do while((tmp2(1) .gt. tol_reinit) .AND.(ii.lt.iter)) ii = ii+1 if(ii.eq.1) then !phi_0 fcc = pcc !phi_0 call mp_level(pcc,ww1,ww2,ww3,ww4) wcase = 1; call weno_phi(pcc,ww1,ww2,wcase) wcase = 2; call weno_phi(pcc,ww3,ww4,wcase) do i=l1s, l1d do j=l1s,l2d sn = pcc(i,j) / sqrt(pcc(i,j)**2.+ww3(i,j)*ddx**2.) sp = max(sn,0.0) ; sm = min(sn,0.0) fcc(i,j) = pcc(i,j) -dt*(sp*ww1(i,j)+sm*ww2(i,j)) enddo enddo !phi_1 == fcc call mp_level(fcc,ww1,ww2,ww3,ww4) wcase = 1; call weno_phi(fcc,ww1,ww2,wcase) wcase = 2; call weno_phi(fcc,ww3,ww4,wcase) do i=l1s, l1d do j=l2s, l2d sn = pcc(i,j) / sqrt(pcc(i,j)**2.+ww3(i,j)*ddx**2.) sp = max(sn,0.0) ; sm = min(sn,0.0) fcc(i,j) = fcc(i,j) -dt*(sp*ww1(i,j)+sm*ww2(i,j)) fcc(i,j) = 3./4.*pcc(i,j)+1./4.*fcc(i,j) enddo enddo !phi_2 == fcc call mp_level(fcc,ww1,ww2,ww3,ww4) wcase = 1; call weno_phi(fcc,ww1,ww2,wcase) wcase = 2; call weno_phi(fcc,ww3,ww4,wcase) do i=l1s, l1d do j=l2s, l2d sn = pcc(i,j) / sqrt(pcc(i,j)**2.+ww3(i,j)*ddx**2.) sp = max(sn,0.0) ; sm = min(sn,0.0) fcc(i,j) = fcc(i,j) -dt*(sp*ww1(i,j)+sm*ww2(i,j)) fcc(i,j) = 1./3.*pcc(i,j)+2./3.*fcc(i,j) enddo enddo !phi_3 == fcc tmp1(1) = 0.0 do i=l1s, l1d do j=l2s, l2d tmp1(1) = tmp1(1)+(fcc(i,j)-pcc(i,j))**2. enddo enddo call globalsum(tmp1,one,tmp2,gcomm) tmp2(1) = sqrt(tmp2(1))/nl1 write(*,*)'[',nloc,'] error in RK3_phi ',tmp2(1) tcc = fcc else call mp_level(tcc,ww1,ww2,ww3,ww4) wcase = 1; call weno_phi(tcc,ww1,ww2,wcase) wcase = 2; call weno_phi(tcc,ww3,ww4,wcase) do i=l1s, l1d do j=l1s,l2d sn = tcc(i,j) / sqrt(tcc(i,j)**2.+ww3(i,j)*ddx**2.) sp = max(sn,0.0) ; sm = min(sn,0.0) fcc(i,j) = tcc(i,j) - dt*(sp*ww1(i,j)+sm*ww2(i,j)) enddo enddo call mp_level(fcc,ww1,ww2,ww3,ww4) wcase = 1; call weno_phi(fcc,ww1,ww2,wcase) wcase = 2; call weno_phi(fcc,ww3,ww4,wcase) do i=l1s, l1d do j=l2s, l2d sn = tcc(i,j) / sqrt(tcc(i,j)**2.+ww3(i,j)*ddx**2.) sp = max(sn,0.0) ; sm = min(sn,0.0) fcc(i,j) = fcc(i,j) - dt*(sp*ww1(i,j)+sm*ww2(i,j)) fcc(i,j) = 3./4.*tcc(i,j)+1./4.*fcc(i,j) enddo enddo call mp_level(fcc,ww1,ww2,ww3,ww4) wcase = 1; call weno_phi(fcc,ww1,ww2,wcase) wcase = 2; call weno_phi(fcc,ww3,ww4,wcase) do i=l1s, l1d do j=l2s, l2d sn = tcc(i,j) / sqrt(tcc(i,j)**2.+ww3(i,j)*ddx**2.) sp = max(sn,0.0) ; sm = min(sn,0.0) fcc(i,j) = fcc(i,j) - dt*(sp*ww1(i,j)+sm*ww2(i,j)) fcc(i,j) = 1./3.*tcc(i,j)+2./3.*fcc(i,j) enddo enddo if(nloc.eq.0) then print *,'in RK3_phi',ii endif tmp1(1) = 0.0 do i=l1s, l1d do j=l2s, l2d tmp1(1) = tmp1(1)+(fcc(i,j)-tcc(i,j))**2. enddo enddo call globalsum(tmp1,one,tmp2,gcomm) tmp2(1) = sqrt(tmp2(1))/nl1 write(*,*)'[',nloc,'] error in RK3_phi ',tmp2(1) tcc = fcc endif enddo return end subroutine RK3_phi subroutine phi_grid(pcc,ccc) use communicators use domain_confg real*8,dimension(l1s:l1d,l2s:l2d) :: ccc real*8,dimension(l1s-3:l1d+3,l2s-3:l2d+3) :: pcc do i=l1s, l1d do j=l2s, l2d ccc(i,j) = pcc(i,j) enddo enddo return end subroutine phi_grid subroutine iphi_grid(ccc,pcc) use communicators use domain_confg real*8,dimension(l1s:l1d,l2s:l2d) :: ccc real*8,dimension(l1s-3:l1d+3,l2s-3:l2d+3) :: pcc do i=l1s, l1d do j=l2s, l2d pcc(i,j) = ccc(i,j) enddo enddo return end subroutine iphi_grid subroutine construct_phi(pcc) use gband_number use communicators use domain_confg use gband_stuff use grid_stuff real*8,dimension(l1s-3:l1d+3,l2s-3:l2d+3) :: pcc real*8 eps, minid,minid1,minid2,minid3,minid4,minid5 real*8 minid6,minid7,minid8 real*8 x1,x2,x3,y1,y2,y3,xx,yy integer nn1,nn2 nn1 = nl1/2 nn2 = nl2*nproc/2 ! eps = 10.**(-8) eps = min(dx,dy) write(*,*)'[',nloc,'] in construct_phi',nn1,nn2 do i=l1s, l1d do j=l2s, l2d tempx1(1:tgb(1)) =abs(float(i-1)*dx-tgbx(1:tgb(1))) tempy1(1:tgb(1)) =abs(float(j-1)*dy-tgby(1:tgb(1))) dd(1:tgb(1)) = sqrt(tempx1**2.+tempy1**2.);minid1=minval(dd) tempx1(1:tgb(1)) =abs(float(i-1)*dx-tgbx(1:tgb(1))) tempy1(1:tgb(1)) =abs(float(j-1)*dy-tgby(1:tgb(1))+nl1*dy) dd(1:tgb(1)) = sqrt(tempx1**2.+tempy1**2.);minid2=minval(dd) tempx1(1:tgb(1)) =abs(float(i-1)*dx-tgbx(1:tgb(1))) tempy1(1:tgb(1)) =abs(float(j-1)*dy-tgby(1:tgb(1))-nl1*dy) dd(1:tgb(1)) = sqrt(tempx1**2.+tempy1**2.);minid3=minval(dd) tempx1(1:tgb(1)) =abs(float(i-1)*dx-tgbx(1:tgb(1))+nl1*dx) tempy1(1:tgb(1)) =abs(float(j-1)*dy-tgby(1:tgb(1))) dd(1:tgb(1)) = sqrt(tempx1**2.+tempy1**2.);minid4=minval(dd) tempx1(1:tgb(1)) =abs(float(i-1)*dx-tgbx(1:tgb(1))-nl1*dx) tempy1(1:tgb(1)) =abs(float(j-1)*dy-tgby(1:tgb(1))) dd(1:tgb(1)) = sqrt(tempx1**2.+tempy1**2.);minid5=minval(dd) minid = min(minid1,minid2,minid3,minid4,minid5) if(pcc(i,j).gt.0) pcc(i,j) = minid if(pcc(i,j).lt.0) pcc(i,j) = -minid if(pcc(i,j).eq.0) pcc(i,j) = 0 enddo enddo return end subroutine construct_phi subroutine weno_grad(pcc,i,j,xy,dir,a) use domain_confg use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff integer i,j,xy,dir real*8, dimension(l1s-3:l1d+3, l2s-3:l2d+3) :: pcc real*8 a,eps,nu1,nu2,nu3,nu4,nu5,s1,s2,s3,a1,a2,a3,w1,w2,w3 eps = 10.**(-6) if(xy.eq.1) then if(dir .gt. 0) then nu1 = (pcc(i-2,j)-pcc(i-3,j))/dx nu2 = (pcc(i-1,j)-pcc(i-2,j))/dx nu3 = (pcc(i, j)-pcc(i-1,j))/dx nu4 = (pcc(i+1,j)-pcc(i, j))/dx nu5 = (pcc(i+2,j)-pcc(i+1,j))/dx endif if(dir .lt. 0) then nu1 = (pcc(i+3,j)-pcc(i+2,j))/dx nu2 = (pcc(i+2,j)-pcc(i+1,j))/dx nu3 = (pcc(i+1,j)-pcc(i, j))/dx nu4 = (pcc(i, j)-pcc(i-1,j))/dx nu5 = (pcc(i-1,j)-pcc(i-2,j))/dx endif s1 = 13./12.*(nu1-2*nu2+ nu3)**2.+ & & 1. / 4.*(nu1-4*nu2+3*nu3)**2. s2 = 13./12.*(nu2-2*nu3+ nu4)**2.+ & & 1. / 4.*(nu2-nu4)**2. s3 = 13./12.*(nu3-2*nu4+ nu5)**2.+ & & 1. / 4.*(3*nu3-4*nu4+nu5)**2. a1 = (1./10.)/(s1+eps)**2. a2 = (6./10.)/(s2+eps)**2. a3 = (3./10.)/(s3+eps)**2. w1 = a1/(a1+a2+a3) w2 = a2/(a1+a2+a3) w3 = a3/(a1+a2+a3) a = w1*( nu1/3.-7*nu2/6.+11.*nu3/6.)+ & & w2*(-nu2/6.+5.*nu3/6.+nu4/3.)+ & & w3*( nu3/3.+5.*nu4/6.-nu5/6.) endif if(xy.eq.2) then if(dir .gt. 0) then nu1 = (pcc(i,j-2)-pcc(i,j-3))/dy nu2 = (pcc(i,j-1)-pcc(i,j-2))/dy nu3 = (pcc(i, j)-pcc(i,j-1))/dy nu4 = (pcc(i,j+1)-pcc(i, j))/dy nu5 = (pcc(i,j+2)-pcc(i,j+1))/dy endif if(dir .le. 0) then nu1 = (pcc(i,j+3)-pcc(i,j+2))/dy nu2 = (pcc(i,j+2)-pcc(i,j+1))/dy nu3 = (pcc(i,j+1)-pcc(i, j))/dy nu4 = (pcc(i, j)-pcc(i,j-1))/dy nu5 = (pcc(i,j-1)-pcc(i,j-2))/dy endif s1 = 13./12.*(nu1-2*nu2+ nu3)**2.+ & & 1. / 4.*(nu1-4*nu2+3*nu3)**2. s2 = 13./12.*(nu2-2*nu3+ nu4)**2.+ & & 1. / 4.*(nu2-nu4)**2. s3 = 13./12.*(nu3-2*nu4+ nu5)**2.+ & & 1. / 4.*(3*nu3-4*nu4+nu5)**2. a1 = (1./10.)/(s1+eps)**2. a2 = (6./10.)/(s2+eps)**2. a3 = (3./10.)/(s3+eps)**2. w1 = a1/(a1+a2+a3) w2 = a2/(a1+a2+a3) w3 = a3/(a1+a2+a3) a = w1*( nu1/3.-7*nu2/6.+11.*nu3/6.)+ & & w2*(-nu2/6.+5.*nu3/6.+nu4/3.)+ & & w3*( nu3/3.+5.*nu4/6.-nu5/6.) endif return end subroutine weno_grad subroutine ynvelocity(cu1,cu2) use domain_confg use physics_stuff use grid_stuff use time_stuff real*8,dimension(l1s:l1d,l2s:l2d) :: cu1,cu2 real*8 f,x0,y0,period, tind,valpha f = 2.*pi/100. x0 = pi ; y0 = pi ; valpha = 0.3 period = 2*pi tind=(time-int(time/period)*period)/period !write(*,*)'[',nloc,'] in ynvelocity ',time,int(time/period),tind do i=l1s, l1d do j=l2s, l2d ! cu1(i,j) = 0.*0.5+ 1.*((j-1)*dy-Ymx/2.) ! cu2(i,j) = 0.*sqrt(3.)/2.+1.*sin((i-1)*dx) ! cu1(i,j) = (50.*f-(j-1)*dy)/10. ! cu2(i,j) = ((i-1)*dx-50.*f)/10. ! cu1(i,j) = 0.-0.*(y-y0)+0.5*sin((j-1.)*dy)-& ! & 0.5*sin((j-1.)*dy)*cos((i-1.)*dx) ! cu2(i,j) = 0.+0.*(x-x0)-0.5*sin((i-1.)*dx)+& ! & 0.5*cos((j-1.)*dy)*sin((i-1.)*dx) ! cu1(i,j) = 2.*pi*cu1(i,j) ! cu2(i,j) = 2.*pi*cu2(i,j) if(tind.lt.(1/2.)) then cu1(i,j) = -2*pi*valpha*sin((j-1)*dy) cu2(i,j) = 0.0 else cu1(i,j) = 0.0 cu2(i,j) = 2*pi*valpha*sin((i-1)*dx) endif cu1(i,j) = ((j-1)*dy-Ymx/2.) cu2(i,j) = sin((i-1)*dx) enddo enddo !print *,'end of ynvelocity' return end subroutine ynvelocity subroutine myvelocity(xx,yy,case,tind,al,& & u,v,ux,uy,vx,vy,uxx,uyy,vxx) real*8 xx,yy,tind,al real*8 u,v,ux,uy,vx,vy,uxx,uyy,vxx integer case real*8 pi pi = 4.*atan(1.0) if(case.eq.1) then !single vortex flow u = pi * sin(yy)*(1-cos(xx)) v =-pi * sin(xx)*(1-cos(yy)) ux = pi*sin(xx)*sin(yy) uy = pi*(1-cos(xx))*cos(yy) vx = -pi*cos(xx)*(1-cos(yy)) vy = -pi*sin(xx)*sin(yy) uxx = pi*cos(xx)*sin(yy) uxy = pi*sin(xx)*cos(yy) uyy = -pi*(1-cos(xx))*sin(yy) vxx = pi*sin(xx)*(1-cos(yy)) endif if(case.eq.2) then !cats eye flow ux = 0.0 ; uy = 1. ; vx = cos(xx) ; vy = 0.0 uxx = 0.0; uxy=0.0 ; uyy=0.0 ; vxx = -sin(xx) endif if(case.eq.3) then !T-dependent K-Flow if(tind.lt.0.5) then u = -2.*pi*al*sin(yy) ; v = 0.0 ux = 0.0 ; uy = -2.*pi*al*cos(yy); vx = 0.0 ; vy = 0.0 uxx= 0.0; uxy =0.0; uyy =2.*pi*al*sin(yy) ; vxx = 0.0 else u = 0.0 ; v = 2.*pi*al*sin(xx) ux = 0.0 ; uy = 0.0; vx = 2.*pi*al*cos(xx); vy = 0.0 uxx= 0.0; uxy =0.0; uyy = 0.0; vxx = -2.*pi*al*sin(xx) endif endif if(case.eq.4) then !rotation flow u = (pi-yy)/10. ; v = (xx-pi)/10. ux = 0.0 ; uy =-1./10. vx = 1./10. ; vy = 0.0 uxx = 0.0; uxy = 0.0; uyy = 0.0 ; vxx = 0.0 endif return end subroutine myvelocity subroutine particledensity(nit) use gband_number use communicators use domain_confg use gband_stuff use grid_stuff use physics_stuff use particle_number include 'mpif.h' integer nit real*8 temp1(1),temp2(1),temp3(1),temp4(1) integer one, root one = 1 root = 0 write(*,*)'[',nloc,'] in pdensity ',mygb(1) temp1(1) = mygb(1) call globalsum(temp1,one,temp2,gcomm) tgb(1) = temp2(1) temp2(1) = totpp(1)+totmp(1) temp3(1) = mypp(2)+mypp(2) call globalsum(temp3,one,temp4,gcomm) write(*,*)'[',nloc,'] total particles',temp4(1),temp2(1) if(nit.eq.1) then pdensity0 = temp4(1)/temp2(1) write(*,*)'[',nloc,'] in p-density', pdensity,pdensity0 endif if(nit.ne.1) then pdensity = temp4(1)/temp2(1) write(*,*)'[',nloc,'] in p-density', pdensity,pdensity0 endif return end subroutine particledensity !************************************************************** subroutine new_surface_tension_force(pcc,ww7,vnx,vny, & & ww1,ww2,ww3,ww4,ww5,ww6) !************************************************************** use domain_phase use steps_stuff use wave_stuff use deriv_stuff use time_stuff use diffuse_stuff use physics_stuff use alias_stuff use grid_stuff use output_stuff use gband_stuff use gband_number use communicators use particle use particle_number include 'mpif.h' real*8, dimension(l1s:l1d,l2s:l2d) :: & & ww1,ww2,ww3,ww4,ww5,ww6,ww7,vnx,vny real*8,dimension(nl1*nl2+nl1*3+nl2*12+36) :: pcc real*8 kx, ky integer idir, itype !ww4 is the curvature call ccurvature1(pcc,ww4) !ww4 is the curvature !ww5 is nx and ww6 is ny call normalvector(pcc,ww5,ww6) !ww5 is nx and ww6 is ny idir = 1 ; itype = 1 call transform(itype,idir,ww7,ww1,ww2) do j=l2s,l2e kx = sk2(j) ! x derivative do i=l1s,l1e-1,2 ky = sk1(i) ! y derivative ww2(i ,j) = -kx * ww1(i+1,j) ww2(i+1,j) = kx * ww1(i, j) ww3(i ,j) = -ky * ww1(i+1,j) ww3(i+1,j) = ky * ww1(i, j) end do end do idir = -1; itype = 1 call transform(itype,idir,ww2,vnx,ww1) call transform(itype,idir,ww3,vny,ww1) do i=l1s, l1e do j=l2s, l2e ww1(i,j) = ww5(i,j) * ww4(i,j) * (vnx(i,j)**2.0d0+vny(i,j)**2.0d0) ww2(i,j) = ww6(i,j) * ww4(i,j) * (vnx(i,j)**2.0d0+vny(i,j)**2.0d0) enddo enddo idir = 1 ; itype = 1 call transform(itype,idir,ww1,ww3,ww6) call transform(itype,idir,ww2,ww4,ww6) do j=l2s,l2e kx = sk2(j) ! x derivative do i=l1s,l1e-1,2 ky = sk1(i) ! y derivative ww5(i ,j) = -ky * ww3(i+1,j) + kx * ww4(i+1,j) ww5(i+1,j) = ky * ww3(i, j) - kx * ww4(i, j) end do end do ! idir = -1 ; itype = 1 ! call transform(itype,idir,ww5,ww6,ww1) ww6 = ww5 do i=l1s, l1e do j=l2s, l2e ww6(i,j) = ww6(i,j)*3.0d0*sqrt(2.0d0*alpha/beta) enddo enddo return end subroutine new_surface_tension_force