/*@@ @file GRHydro_BoundariesM.F90 @date Aug 30, 2010 @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_InitSymBoundM @date Aug 30, 2010 @author Joshua Faber, 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_InitSymBoundM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer, dimension(3) :: sym integer :: ierr 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(clean_divergence.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::psidc") endif !!$ handle multiple tracer variables if(evolve_tracer.eq.1) then call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_tracers") call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_cons_tracers") endif if(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then call SetCartSymGN(ierr, cctkGH, sym, "HydroBase::Y_e") call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::Y_e_con") 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]") 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]") 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]") call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]") end subroutine GRHydro_InitSymBoundM /*@@ @routine GRHydro_BoundariesM @date Aug 30, 2010 @author @desc Calls the appropriate boundary routines @enddesc @calls @calledby @history @endhistory @@*/ subroutine GRHydro_BoundariesM(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") 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 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(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) 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 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") 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 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(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) 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 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_BoundariesM /*@@ @routine GRHydro_OutflowBoundariesM @date Aug 30, 2010 @author Joshua Faber, 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_OutflowBoundariesM(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) 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 psi4pt = 1.d0 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)) 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)) 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) 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 psi4pt = 1.d0 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)) 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)) end if end do end do end do end if end if end subroutine GRHydro_OutflowBoundariesM