subroutine hydro_3d (mype)
!-------------------------------------------------------------------
!
! Step 1
!
!-------
! Include the amr package's header files.
use paramesh_dimensions
use physicaldata
use tree
use workspaceuse paramesh_interfaces
use paramesh_mpi_interfaces
include 'mpif.h' ! Include header files for user's code. include 'user_header_file.fh' ! Declare any local variables required. real dthalf real rho(1:nxb+nguard,1:nzb+nguard,1:nzb+nguard) real vx(1:nxb+nguard,1:nzb+nguard,1:nzb+nguard) real vy(1:nxb+nguard,1:nzb+nguard,1:nzb+nguard) real vz(1:nxb+nguard,1:nzb+nguard,1:nzb+nguard) real ener(1:nxb+nguard,1:nzb+nguard,1:nzb+nguard) real xfluxes(5,1:nxb+nguard+1,1:nzb+nguard,1:nzb+nguard) real yfluxes(5,1:nxb+nguard,1:nzb+nguard+1,1:nzb+nguard) real zfluxes(5,1:nxb+nguard,1:nzb+nguard,1:nzb+nguard+1) real src(5,1:nxb+nguard,1:nzb+nguard,1:nzb+nguard) integer mype,iopt,nlayers,jface,lb !------------------------------------------------------------------- ! ! Step 2 ! !------- ! Make a temporary copy of the solution at the beginning of the timestep. t_unk(:,:,:,:,:) = unk(:,:,:,:,:) ! Integrate forward through a half timestep dthalf = dt*.5 ! loop over leaf blocks on this processor if(lnblocks.gt.0) then do lb=1,lnblocks if(nodetype(lb).eq.1) then ! copy solution data for this block from amr datastructure to user's local ! variables rho(:,:,:) = unk(1,:,:,:,lb) vx(:,:,:) = unk(2,:,:,:,lb) vy(:,:,:) = unk(3,:,:,:,lb) vz(:,:,:) = unk(4,:,:,:,lb) ener(:,:,:) = unk(5,:,:,:,lb) ! Compute required source terms on this block call source(src,rho,vx,vy,vz,ener) ! Advance the solution on this block through dt/2 call advance(dthalf,rho,vx,vy,vz,ener, . xfluxes,yfluxes,zfluxes,src) endif enddo endif !------------------------------------------------------------------- ! ! Step 3 ! !------- ! Fill the guard cells. iopt = 1 nlayers = nguard call amr_guardcell(mype,iopt,nlayers) !------------------------------------------------------------------- ! ! Step 4 ! !------- ! Integrate forward through a full timestep ! loop over leaf blocks on this processor if(lnblocks.gt.0) then do lb=1,lnblocks if(nodetype(lb).eq.1) then ! copy solution data for this block from amr datastructure to user's local ! variables rho(:,:,:) = unk(1,:,:,:,lb) vx(:,:,:) = unk(2,:,:,:,lb) vy(:,:,:) = unk(3,:,:,:,lb) vz(:,:,:) = unk(4,:,:,:,lb) ener(:,:,:) = unk(5,:,:,:,lb) ! Compute required source terms on this block call source(src,rho,vx,vy,vz,ener) ! Restore original solution unk(1,:,:,:,lb) = t_unk(1,:,:,:,lb) unk(2,:,:,:,lb) = t_unk(2,:,:,:,lb) unk(3,:,:,:,lb) = t_unk(3,:,:,:,lb) unk(4,:,:,:,lb) = t_unk(4,:,:,:,lb) unk(5,:,:,:,lb) = t_unk(5,:,:,:,lb) ! Advance the solution on this block through dt/2 call advance(dt,rho,vx,vy,vz,ener, . xfluxes,yfluxes,zfluxes,src) ! Capture fluxes at block boundaries flux_x(:,1,:,:,lb) = xfluxes(:,nguard,:,:) flux_x(:,2,:,:,lb) = xfluxes(:,nxb+nguard+1,:,:) flux_y(:,:,1,:,lb) = yfluxes(:,:,nguard,:) flux_y(:,:,2,:,lb) = yfluxes(:,:,nyb+nguard+1,:) flux_z(:,:,:,1,lb) = zfluxes(:,:,:,nguard) flux_z(:,:,:,2,lb) = zfluxes(:,:,:,nzb+nguard+1) endif enddo endif !------------------------------------------------------------------- ! ! Step 5 ! !------- ! Modify the boundary fluxes ! ! note: the first thing that happens in amr_flux_conserve is that the ! arrays flux_x, flux_y and flux_z are copied into temporary storage ! arrays tflux_x, tflux_y and tflux_z respectively. nsub = 1 call amr_flux_conserve(mype,nsub) ! loop over leaf blocks on this processor if(lnblocks.gt.0) then do lb=1,lnblocks if(nodetype(lb).eq.1) then dx = bsize(1,lb)/real(nxb) dy = bsize(2,lb)/real(nyb) dz = bsize(3,lb)/real(nzb) area_yz = dy*dz area_xz = dx*dz area_xy = dx*dy vol = dx*dy*dz rvol = 1./vol ! Update the solution immediately inside the block boundaries to reflect this ! modification of the fluxes. Note it is assumed here that the arrays flux_x, ! tflux_x, etc are fluxes per unit area. unk(:,1+nguard,:,:,lb) = unk(:,1+nguard,:,:,lb) + . (flux_x(:,1,:,:,lb) - tflux_x(:,1,:,:,lb))*area_yz*rvol*dt unk(:,nxb+nguard,:,:,lb) = unk(:,nxb+nguard,:,:,lb) + . (tflux_x(:,2,:,:,lb) - flux_x(:,2,:,:,lb))*area_yz*rvol*dt unk(:,:,1+nguard,:,lb) = unk(:,:,1+nguard,:,lb) + . (flux_y(:,:,1,:,lb) - tflux_y(:,:,1,:,lb))*area_xz*rvol*dt unk(:,nyb+nguard,:,lb) = unk(:,:,nyb+nguard,:,lb) + . (tflux_y(:,:,2,:,lb) - flux_y(:,:,2,:,lb))*area_xz*rvol*dt unk(:,:,:,1+nguard,lb) = unk(:,:,:,1+nguard,lb) + . (flux_z(:,:,:,1,lb) - tflux_z(:,:,:,1,lb))*area_xy*rvol*dt unk(:,:,:,nzb+nguard,lb) = unk(:,:,:,nzb+nguard,lb) + . (tflux_z(:,:,:,2,lb) - flux_z(:,:,:,2,lb))*area_xy*rvol*dt endif enddo endif !------------------------------------------------------------------- return end program hydro_3d ! User supplied source routine. subroutine source(src,rho,vx,vy,vz,ener) ! ! Arguments : ! src real array source terms ! rho real array density ! vx real array x-component of velocity ! vy real array y-component of velocity ! vz real array z-component of velocity ! ener real array energy density ! Include header files for user's code. include 'user_header_file.fh' ! defines parameters nx,ny,nz real rho(1:nx+4,1:nz+4,1:nz+4) real vx (1:nx+4,1:nz+4,1:nz+4) real vy (1:nx+4,1:nz+4,1:nz+4) real vz (1:nx+4,1:nz+4,1:nz+4) real ener(1:nx+4,1:nz+4,1:nz+4) real src(5,1:nx+4,1:nz+4,1:nz+4) do k=3,nz+2 do j=3,ny+2 do i=3,nx+2 src(1,i,j,k) = ....... src(2,i,j,k) = ....... src(3,i,j,k) = ....... src(4,i,j,k) = ....... src(5,i,j,k) = ....... enddo enddo enddo return end subroutine source ! User supplied solution advance routine. subroutine advance(dt,rho,vx,vy,vz,ener,xfluxes,yfluxes,zfluxes,src) ! This routine assumes 2 guard cells at each boundary. ! ! Arguments : ! dt real timestep ! rho real array density ! vx real array x-component of velocity ! vy real array y-component of velocity ! vz real array z-component of velocity ! ener real array energy density ! xfluxes real array flux densities multiplied by area of ! cell face across x faces ! yfluxes real array flux densities multiplied by area of ! cell face across y faces ! zfluxes real array flux densities multiplied by area of ! cell face across z faces ! src real array source terms ! Include header files for user's code. include 'user_header_file.fh' ! defines parameters nx,ny,nz real rho(1:nx+4,1:nz+4,1:nz+4) real vx (1:nx+4,1:nz+4,1:nz+4) real vy (1:nx+4,1:nz+4,1:nz+4) real vz (1:nx+4,1:nz+4,1:nz+4) real ener(1:nx+4,1:nz+4,1:nz+4) real xfluxes(5,1:nx+5,1:nz+4,1:nz+4) real yfluxes(5,1:nx+4,1:nz+5,1:nz+4) real zfluxes(5,1:nx+4,1:nz+4,1:nz+5) real src(5,1:nx+4,1:nz+4,1:nz+4) real dt ! Compute fluxes across x-faces of grid cells do k=3,nz+2 do j=3,ny+2 do i=3,nx+2+1 xfluxes(1,i,j,k) = .... xfluxes(2,i,j,k) = .... xfluxes(3,i,j,k) = .... xfluxes(4,i,j,k) = .... xfluxes(5,i,j,k) = .... enddo enddo enddo ! Compute fluxes across y-faces of grid cells do k=3,nz+2 do j=3,ny+2+1 do i=3,nx+2 yfluxes(1,i,j,k) = .... yfluxes(2,i,j,k) = .... yfluxes(3,i,j,k) = .... yfluxes(4,i,j,k) = .... yfluxes(5,i,j,k) = .... enddo enddo enddo ! Compute fluxes across z-faces of grid cells do k=3,nz+2+1 do j=3,ny+2 do i=3,nx+2 zfluxes(1,i,j,k) = .... zfluxes(2,i,j,k) = .... zfluxes(3,i,j,k) = .... zfluxes(4,i,j,k) = .... zfluxes(5,i,j,k) = .... enddo enddo enddo ! Update solution using these fluxes do k=3,nz+2 do j=3,ny+2 do i=3,nx+2 dx = bsize(1,lb)/real(nxb) dy = bsize(2,lb)/real(nyb) dz = bsize(3,lb)/real(nzb) area_yz = dy*dz area_xz = dx*dz area_xy = dx*dy vol = dx*dy*dz rvol = 1./vol rho(i,j,k) = rho(i,j,k) + rvol*dt*( . (xfluxes(1,i,j,k) - xfluxes(1,i+1,j,k))*area_yz + . (yfluxes(1,i,j,k) - yfluxes(1,i,j+1,k))*area_zx + . (zfluxes(1,i,j,k) - zfluxes(1,i,j,k+1))*area_xy ) . + src(1,i,j,k)*dt vx(i,j,k) = vx(i,j,k) + rvol*dt*( . (xfluxes(2,i,j,k) - xfluxes(2,i+1,j,k))*area_yz + . (yfluxes(2,i,j,k) - yfluxes(2,i,j+1,k))*area_zx + . (zfluxes(2,i,j,k) - zfluxes(2,i,j,k+1))*area_xy ) . + src(2,i,j,k)*dt vy(i,j,k) = vy(i,j,k) + rvol*dt*( . (xfluxes(3,i,j,k) - xfluxes(3,i+1,j,k))*area_yz + . (yfluxes(3,i,j,k) - yfluxes(3,i,j+1,k))*area_zx + . (zfluxes(3,i,j,k) - zfluxes(3,i,j,k+1))*area_xy ) . + src(3,i,j,k)*dt vz(i,j,k) = vz(i,j,k) + rvol*dt*( . (xfluxes(4,i,j,k) - xfluxes(4,i+1,j,k))*area_yz + . (yfluxes(4,i,j,k) - yfluxes(4,i,j+1,k))*area_zx + . (zfluxes(4,i,j,k) - zfluxes(4,i,j,k+1))*area_xy ) . + src(4,i,j,k)*dt ener(i,j,k) = ener(i,j,k) + rvol*dt*( . (xfluxes(5,i,j,k) - xfluxes(5,i+1,j,k))*area_yz + . (yfluxes(5,i,j,k) - yfluxes(5,i,j+1,k))*area_zx + . (zfluxes(5,i,j,k) - zfluxes(5,i,j,k+1))*area_xy ) . + src(5,i,j,k)*dt enddo enddo enddo return end subroutine advance