subroutine advect_muscl(mype,dt,time,nprocs,loop)
!---------------------------------------
!
! This routine assumes that NO_PERMANENT_GUARDCELLS is
! undefined in paramesh_preprocessor.fh.
! Include files for amr package
use paramesh_dimensions
use physicaldata
use tree
use paramesh_interfaces
use paramesh_mpi_interfaces
include 'pointers.fh'
!---------------------------------------
implicit none
include 'mpif.h'
real :: dtp
real :: rho_s1(maxblocks),rho_s2(maxblocks)
real :: rho_sum1,rho_sum2
save :: dtp,rho_sum1,rho_sum2
integer :: mype,loop,iopt,nlayers,nprocs,nsub
logical :: l_bc
!-------------------------------------------------------------------------------
! USERS CODE BLOCK
integer, parameter :: mx = nxb+2*nguard !<<< EDITED
integer, parameter :: mxp = mx+1
integer, parameter :: mxm = mx-1
! number of primitives and unknowns
integer, parameter :: neq = nvar !<<< EDITED
integer :: i,j
! Logical Variables
logical :: lstp ! if true use steepner
! Global Variables
real, dimension(neq,mx) :: vvg ! primitives
real, dimension(neq,mx) :: dvg ! incremental change to primitives
real, dimension(neq,mxp):: flx ! grid cell interface fluxes
real, dimension(mx ) :: cs ! sound speed
real, dimension(neq,mx) :: wp, w , wm
real, dimension(neq,mx) :: dw, dwdum
real, dimension(neq,mx) :: dvl,dvl1
real, dimension(neq,mx) :: dvr,dvr1,dvgd
real, dimension(neq,mx) :: eigv
real, dimension(mx) :: evmax, evmin,dtdx
real, dimension(neq) :: chi
real, dimension(neq,neq,mx) :: lhev,rhev
real, dimension(neq,mx) :: scr1 ,scr2 ,scr3 ,scr4
real, dimension(neq,mx) :: scr5 ,scr6 ,scr7 ,scr8
real, dimension(neq,mx) :: scr9 ,scr10,scr11,scr12
real, dimension(neq,mx) :: scr13,scr14,scr15,scr16
common/minima/rmin,pmin,smin
common/control/crn,fact,nriem,ibc
real :: time, crn, fact, pmin, rmin, smin, dt, dth
integer :: ibc
integer :: nriem
integer :: nx,nxp
! Boundary Condition Variables
integer, dimension(2) :: lbloc,bctype
! Miscellaneous Local Variables
real :: dtl
real :: xmin,xmax,delx,delxi
integer :: lblock
save dtl
! Mesh Variables
real, dimension(mx) :: rvol ! inverse control volume at i,j
real, dimension(mx) :: xc ! x-position of control volume center
real, dimension(mxp) :: xm ! xm(i) is the cell interface position
! to the left of the cell center
real, dimension(mx) :: dxm
real, dimension(mx) :: rdx
real, dimension(mx) :: dxc ! dx is the distance from xm(i) to
! xm(i+1)
real, dimension(mx) :: rdxc ! rdxc is the inverse of dxc
chi(1)=2.
chi(2)=1.
chi(3)=2.
nx = mx
nxp = nx+1
lstp = .true.
! END OF USERS CODE BLOCK
!-------------------------------------------------------------------------------
rho_s1(:) = 0.
rho_s2(:) = 0.
rho_sum1 = 0.
!------------------------------------
! Section 1. ! Begin Global Timestep Determination !------------------------------------ dtp = 1.e30 if(lnblocks.gt.0) then do lblock=1,lnblocks if(nodetype(lblock).eq.1) then ! set up pointers uug => unk(:,:,1,1,lblock) xmax = bnd_box(2,1,lblock) xmin = bnd_box(1,1,lblock) delx = (xmax - xmin)/real(nxb) !------------------------------------------------------------------------------- ! USERS CODE BLOCK delxi = 1./delx do i=1,nx rdx(i) = delxi enddo ! determine primitives from conserved values call eos(uug,vvg,cs,pmin,rmin,mx,nx,neq) call dtset(vvg,cs,mx,nx,neq,crn,rdx,dtl) ! END OF USERS CODE BLOCK !------------------------------------------------------------------------------- ! Keep track of the minimum timestep amongst grid blocks on this processor dtp = min(dtp,dtl) do i=1+nguard,nxb+nguard rho_s1(lblock) = rho_s1(lblock) . +unk(1,i,1,1,lblock)*delx enddo rho_sum1 = rho_sum1 + rho_s1(lblock) endif enddo endif call comm_real_sum_to_all(rho_sum1,rho_sum1) ! find smallest hydro timestep across all processors call comm_real_min_to_all(dtp,dtp) dt = dtp ! write(*,*) 'Timestep = ',dt !--------------------------------------- ! Global Timestep Now Determined. ! End of Section 1. !--------------------------------------- !--------------------------------------- ! Section 2. ! Now compute fluxes using Muscl scheme. !--------------------------------------- rho_sum2 = 0. if(lnblocks.gt.0) then do lblock=1,lnblocks if(nodetype(lblock).eq.1) then ! set up pointers uug => unk(:,:,1,1,lblock) xmax = bnd_box(2,1,lblock) xmin = bnd_box(1,1,lblock) delx = (xmax - xmin)/real(nxb) !------------------------------------------------------------------------------- ! USERS CODE BLOCK call eos(uug,vvg,cs,pmin,rmin,mx,nx,neq) call ngrid ( xmin,xmax,delx,rdx, . xc ,xm ,dxc ,dxm, . rdxc,fact,rvol, . mx ,mxp ,nx , nxp ) dtdx(:nx) = dt*rdxc(:nx) call cfd_1d( . uug ,vvg ,dvg ,flx ,wp ,w , . wm ,dw ,dwdum,dxc ,dvl ,dvr , . dvl1 ,dvr1 ,dvgd ,evmax,evmin , . lbloc,bctype, . eigv,chi ,cs ,xc ,rvol , . scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8, . scr9,scr10,scr11,scr12,scr13,scr14,scr15, . lhev ,rhev , . lstp ,time ,dtdx,dt , . rmin ,pmin ,smin ,mx ,nx ,neq ,nriem, . ibc) ! END OF USERS CODE BLOCK !------------------------------------------------------------------------------- ! Capture fluxes applied to solution at block boundaries ! First left x block boundary. Since this is a 1D code the appropriate ! indeces for flux_x are ! index 1 : meaning apply to all cell centered variables ! index 2 1 indicates left x block boundary ! index 3 1 since in 1D there is only one cell in the y direction ! index 4 1 since in 1D there is only one cell in the z direction ! index 5 lblock the local block number flux_x(:,1,1,1,lblock) = flx(:,1+nguard) ! Now right x block boundary flux_x(:,2,1,1,lblock) = flx(:,1+nxb+nguard) endif enddo endif !--------------------------------------- ! End of Section 2. !--------------------------------------- !---------------------------------------------------------------- ! Now patch up fluxes at block boundaries to ensure conservation. !---------------------------------------------------------------- nsub = 1 call amr_flux_conserve(mype,nsub)
!---------------------------------------------------------------- ! Section 3. ! Final step, apply patched up fluxes to advance solution. !---------------------------------------------------------------- if(lnblocks.gt.0) then do lblock=1,lnblocks if(nodetype(lblock).eq.1) then ! set up pointers uug => unk(:,:,1,1,lblock) xmax = bnd_box(2,1,lblock) xmin = bnd_box(1,1,lblock) delx = (xmax - xmin)/real(nxb) ! Is this block at an external boundary? l_bc = .false. if(neigh(1,1,lblock).le.-20.or. . neigh(1,2,lblock).le.-20) l_bc=.true. !------------------------------------------------------------------------------- ! USERS CODE BLOCK call ngrid ( xmin,xmax,delx,rdx, . xc ,xm ,dxc ,dxm, . rdxc,fact,rvol, . mx ,mxp ,nx , nxp ) ! Correct the first interior cells at each block boundary using the adjusted ! fluxes returned from amr_flux_conserve i = 1+nguard uug(:neq,i)= uug(:neq,i) . -( tflux_x(:neq,1,1,1,lblock) . - flux_x(:neq,1,1,1,lblock) )*rvol(i) i= nxb+nguard uug(:neq,i)= uug(:neq,i) . -( flux_x(:neq,2,1,1,lblock) . -tflux_x(:neq,2,1,1,lblock) )*rvol(i) if(l_bc) call bcset(uug,mx,nx,neq,rmin,pmin,ibc) !<<< EDITED ! END OF USERS CODE BLOCK !------------------------------------------------------------------------------- endif enddo endif !--------------------------------------- ! End of Section 3. !--------------------------------------- if(lnblocks.gt.0) then do lblock=1,lnblocks if(nodetype(lblock).eq.1) then xmax = bnd_box(2,1,lblock) xmin = bnd_box(1,1,lblock) delx = (xmax - xmin)/real(nxb) do i=1+nguard,nxb+nguard rho_s2(lblock) = rho_s2(lblock) . +unk(1,i,1,1,lblock)*delx enddo rho_sum2 = rho_sum2 + rho_s2(lblock) endif enddo endif time = time + dt iopt=1 nlayers = nguard call amr_guardcell(mype,iopt,nlayers) call comm_real_sum_to_all(rho_sum2,rho_sum2) ! call MPI_BARRIER(MPI_COMM_WORLD, ierr) ! if(mype.eq.0) then ! write(*,*) ' ' ! write(*,*) 'Conservation Check' ! write(*,*) 'rho_sum before dt = ',rho_sum1 ! write(*,*) 'rho_sum after dt = ',rho_sum2 ! write(*,*) ' ' ! endif
return end subroutine advect_muscl