/*@@ @file GRHydro_Boundaries.F90 @date Sat Jan 26 01:01:14 2002 @author @desc The two routines for dealing with boundary conditions. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" #include "GRHydro_Macros.h" #include "util_Table.h" #define velx(i,j,k) vel(i,j,k,1) #define vely(i,j,k) vel(i,j,k,2) #define velz(i,j,k) vel(i,j,k,3) #define sx(i,j,k) scon(i,j,k,1) #define sy(i,j,k) scon(i,j,k,2) #define sz(i,j,k) scon(i,j,k,3) #define Bvecx(i,j,k) Bvec(i,j,k,1) #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) /*@@ @routine GRHydro_InitSymBound @date Sat Jan 26 01:03:04 2002 @author Ian Hawke @desc Sets up the symmetries at the boundaries of the hydrodynamical variables. @enddesc @calls @calledby @history Direct translation of routines from GR3D, GRAstro_Hydro, written by Mark Miller, or WaveToy routines, or... @endhistory @@*/ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer, dimension(3) :: sym integer :: ierr = 0 integer :: itracer character(len=100) tracername character(len=100) tracerindex sym(1) = 1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::rho") call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::press") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::dens") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::tau") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::w_lorentz") call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::eps") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::GRHydro_C2P_failed") if(evolve_mhd.ne.0.and.clean_divergence.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::psidc") endif !!$ handle multiple tracer variables if(evolve_tracer.ne.0) then call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_tracers") call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_cons_tracers") endif if(evolve_y_e.ne.0) then call SetCartSymGN(ierr, cctkGH, sym, "HydroBase::Y_e") call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::Y_e_con") endif if(evolve_temper.ne.0) then call SetCartSymGN(ierr, cctkGH, sym, "HydroBase::temperature") call SetCartSymGN(ierr, cctkGH, sym, "HydroBase::entropy") endif sym(1) = -1 sym(2) = 1 sym(3) = 1 call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]") if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]") sym(1) = 1 sym(2) = -1 sym(3) = 1 call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]") if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") sym(1) = 1 sym(2) = 1 sym(3) = -1 call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]") if(evolve_mhd.ne.0)call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]") end subroutine GRHydro_InitSymBound /*@@ @routine GRHydro_Boundaries @date Sat Jan 26 01:04:04 2002 @author @desc Calls the appropriate boundary routines @enddesc @calls @calledby @history @endhistory @@*/ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer, dimension(3) :: sw integer :: ierr CCTK_REAL :: pi = 3.141569d0 integer :: i,j,k CCTK_INT, parameter :: faces=CCTK_ALL_FACES CCTK_INT, parameter :: ione=1 sw = GRHydro_stencil !!$Flat boundaries if required if (CCTK_EQUALS(bound,"flat")) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::dens", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::tau", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::scon", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::w_lorentz", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::rho", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::press", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::eps", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::vel", "Flat") if(evolve_mhd.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Bvec", "Flat") if(clean_divergence.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::psidc", "Flat") endif endif if(evolve_tracer.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::GRHydro_tracers", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::GRHydro_cons_tracers", "Flat") endif if(evolve_y_e.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Y_e", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Y_e_con", "Flat") endif if(evolve_temper.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::temperature", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::entropy", "Flat") endif endif if (CCTK_EQUALS(bound,"none")) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::dens", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::tau", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::scon", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::w_lorentz", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::rho", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::press", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::eps", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::vel", "None") if(evolve_mhd.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Bvec", "None") if(clean_divergence.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::psidc", "None") endif endif if(evolve_tracer.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::GRHydro_tracers", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::GRHydro_cons_tracers", "None") endif if(evolve_y_e.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::Y_e", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Y_e_con", "None") endif if(evolve_temper.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::temperature", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::entropy", "None") endif end if if (CCTK_EQUALS(bound,"scalar")) then call CCTK_WARN(0, "Until somebody uses this I see no reason to support it") end if if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") end subroutine GRHydro_Boundaries /*@@ @routine GRHydro_OutflowBoundaries @date Tue May 17 16:58:06 2005 @author Ian Hawke @desc Set outflow boundaries over only part of the domain. This is designed to be used with GZPatchSystem. @enddesc @calls @calledby @history @endhistory @@*/ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer, dimension(3) :: sw integer :: ierr integer :: i,j,k CCTK_REAL, dimension(3) :: posn CCTK_REAL :: det,psi4pt sw = GRHydro_stencil if (r(1,1,1) < r(1,1,cctk_lsh(3))) then if (cctk_bbox(6) .ne. 0) then do k = cctk_lsh(3) - sw(3), cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) posn(1) = x(i,j,k) posn(2) = y(i,j,k) posn(3) = z(i,j,k) if (dot_product(outflowboundary_normal, posn) > 0.d0) then rho(i,j,k) = rho(i,j,k-1) velx(i,j,k) = velx(i,j,k-1) vely(i,j,k) = vely(i,j,k-1) velz(i,j,k) = velz(i,j,k-1) eps(i,j,k) = eps(i,j,k-1) press(i,j,k) = press(i,j,k-1) w_lorentz(i,j,k) = w_lorentz(i,j,k-1) if(evolve_mhd.ne.0) then Bvecx(i,j,k) = Bvecx(i,j,k-1) Bvecy(i,j,k) = Bvecy(i,j,k-1) Bvecz(i,j,k) = Bvecz(i,j,k-1) if(clean_divergence.ne.0) then psidc(i,j,k) = psidc(i,j,k-1) endif endif psi4pt = 1.0d0 det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) if(evolve_mhd.ne.0) then call prim2conM(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) else call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) endif end if end do end do end do end if else if (cctk_bbox(5) .ne. 0) then do k = sw(3), 1, -1 do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) posn(1) = x(i,j,k) posn(2) = y(i,j,k) posn(3) = z(i,j,k) if (dot_product(outflowboundary_normal, posn) > 0.d0) then rho(i,j,k) = rho(i,j,k+1) velx(i,j,k) = velx(i,j,k+1) vely(i,j,k) = vely(i,j,k+1) velz(i,j,k) = velz(i,j,k+1) eps(i,j,k) = eps(i,j,k+1) press(i,j,k) = press(i,j,k+1) w_lorentz(i,j,k) = w_lorentz(i,j,k+1) if(evolve_mhd.ne.0) then Bvecx(i,j,k) = Bvecx(i,j,k+1) Bvecy(i,j,k) = Bvecy(i,j,k+1) Bvecz(i,j,k) = Bvecz(i,j,k+1) if(clean_divergence.ne.0) then psidc(i,j,k) = psidc(i,j,k+1) endif endif psi4pt = 1.0d0 det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) if(evolve_mhd.ne.0) then call prim2conM(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),& rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) else call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) endif end if end do end do end do end if end if end subroutine GRHydro_OutflowBoundaries