From 74fb1e6ea34d6e03a35ff6c158f455c39904bf5a Mon Sep 17 00:00:00 2001 From: bmundim Date: Sun, 2 May 2010 20:59:32 +0000 Subject: file/parameter string replacement from whisky to GRHydro git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@112 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Boundaries.F90 | 326 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 326 insertions(+) create mode 100644 src/GRHydro_Boundaries.F90 (limited to 'src/GRHydro_Boundaries.F90') diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 new file mode 100644 index 0000000..5e8c093 --- /dev/null +++ b/src/GRHydro_Boundaries.F90 @@ -0,0 +1,326 @@ + /*@@ + @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 "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) + + /*@@ + @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 + + 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") + +!!$ 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 + + 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]") + + 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]") + + 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]") + +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_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 + + 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_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 + + 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 (conformal_state .eq. 0) then + psi4pt = 1d0 + call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) + else + psi4pt = psi(i,j,k)**4 + call SpatialDeterminant(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) + end if + 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)) + + 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 (conformal_state .eq. 0) then + psi4pt = 1d0 + call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) + else + psi4pt = psi(i,j,k)**4 + call SpatialDeterminant(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) + end if + 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)) + + end if + + end do + end do + end do + + end if + + end if + +end subroutine GRHydro_OutflowBoundaries + + + -- cgit v1.2.3