From e0dc2af4862d5ddb874328bd097f7f516231dd8c Mon Sep 17 00:00:00 2001 From: rhaas Date: Thu, 15 Sep 2011 16:49:53 +0000 Subject: add Multipatch support to GRHydro * not all features of GRHydro are supported yet, in particular only the HLLE solver supports Mulitpatch yet. Original commit by Christian Reisswig and Christian Ott git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@273 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Boundaries.F90 | 24 +- src/GRHydro_Con2Prim.F90 | 114 ++++----- src/GRHydro_ENOReconstruct_drv.F90 | 12 +- src/GRHydro_EoSChangeGamma.F90 | 22 +- src/GRHydro_HLLE.F90 | 48 ++-- src/GRHydro_Jacobian_state.c | 32 +++ src/GRHydro_PPMMReconstruct_drv.F90 | 24 +- src/GRHydro_PPMReconstruct_drv.F90 | 18 +- src/GRHydro_Prim2Con.F90 | 224 +++++++++--------- src/GRHydro_RegisterVars.cc | 8 +- src/GRHydro_Source.F90 | 432 +++++++++++++++++------------------ src/GRHydro_TVDReconstruct_drv.F90 | 24 +- src/GRHydro_TransformTensorBasis.F90 | 211 +++++++++++++++++ src/GRHydro_UpdateMask.F90 | 76 +++--- src/make.code.defn | 13 +- 15 files changed, 766 insertions(+), 516 deletions(-) create mode 100644 src/GRHydro_Jacobian_state.c create mode 100644 src/GRHydro_TransformTensorBasis.F90 (limited to 'src') diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 index 5bfd457..0a1b3d8 100644 --- a/src/GRHydro_Boundaries.F90 +++ b/src/GRHydro_Boundaries.F90 @@ -15,15 +15,15 @@ #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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) +#define Bvecx(i,j,k) lBvec(i,j,k,1) +#define Bvecy(i,j,k) lBvec(i,j,k,2) +#define Bvecz(i,j,k) lBvec(i,j,k,3) #define Bconsx(i,j,k) Bcons(i,j,k,1) #define Bconsy(i,j,k) Bcons(i,j,k,2) #define Bconsz(i,j,k) Bcons(i,j,k,3) @@ -94,7 +94,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(2) = 1 sym(3) = 1 - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[0]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]") if(evolve_mhd.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]") @@ -105,7 +105,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(2) = -1 sym(3) = 1 - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[1]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]") if(evolve_mhd.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") @@ -116,7 +116,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(2) = 1 sym(3) = -1 - call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]") + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[2]") call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]") if(evolve_mhd.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]") @@ -182,10 +182,10 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::eps", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "Flat") + "GRHydro::lvel", "Flat") if(evolve_mhd.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") + "GRHydro::lBvec", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Bcons", "Flat") if(clean_divergence.ne.0) then @@ -237,7 +237,7 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) "HydroBase::vel", "None") if(evolve_mhd.ne.0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") + "GRHydro::lBvec", "None") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Bcons", "None") if(clean_divergence.ne.0) then diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index d3c9732..3c51302 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -101,11 +101,11 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) epsnegative = .false. - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) + gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),& + gbc(i,j,k),gcc(i,j,k)) !!$ if (det < 0.e0) then !!$ call CCTK_WARN(0, "nan produced (det negative)") @@ -135,7 +135,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 - vel(i,j,k,:) = 0.d0 + lvel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 if(evolve_temper.ne.0) then @@ -148,7 +148,6 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& press(i,j,k),keyerr,anyerr) else - ! w_lorentz=1, so the expression for tau reduces to: keytemp = 0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) @@ -157,6 +156,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) endif + ! w_lorentz=1, so the expression for tau reduces to: tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k)) - dens(i,j,k) cycle @@ -165,15 +165,15 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) if(evolve_temper.eq.0) then call Con2Prim_pt(GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2), & - scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), & - vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & + scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), & + lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) else call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& - scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vel(i,j,k,1),& - vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),lvel(i,j,k,1),& + lvel(i,j,k,2),lvel(i,j,k,3),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & @@ -186,7 +186,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) local_perc_ptol = GRHydro_perc_ptol*100.0d0 call Con2Prim_pt_hot(i,j,k,GRHydro_eos_handle, dens(i,j,k),scon(i,j,k,1),& scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),& - vel(i,j,k,1),vel(i,j,k,2), vel(i,j,k,3),eps(i,j,k),& + lvel(i,j,k,1),lvel(i,j,k,2), lvel(i,j,k,3),eps(i,j,k),& temperature(i,j,k),y_e(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & @@ -237,7 +237,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) rho(i,j,k) = GRHydro_rho_min scon(i,j,k,:) = 0.d0 - vel(i,j,k,:) = 0.d0 + lvel(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& @@ -259,8 +259,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) !$OMP END CRITICAL call Con2Prim_ptPolytype(GRHydro_polytrope_handle, & dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & - tau(i,j,k), rho(i,j,k), vel(i,j,k,1), vel(i,j,k,2), & - vel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & + tau(i,j,k), rho(i,j,k), lvel(i,j,k,1), lvel(i,j,k,2), & + lvel(i,j,k,3), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k), & uxx, uxy, uxz, uyy, uyz, uzz, det, x(i,j,k), y(i,j,k), & z(i,j,k), r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) #endif @@ -1048,18 +1048,18 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) epsnegative = .false. @@ -1251,11 +1251,11 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) + gaa(i,j,k),gab(i,j,k),gac(i,j,k),gbb(i,j,k),& + gbc(i,j,k),gcc(i,j,k)) if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -1277,8 +1277,8 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) call Con2Prim_ptPolytype(GRHydro_eos_handle, dens(i,j,k),& scon(i,j,k,1),scon(i,j,k,2), & - scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vel(i,j,k,1),vel(i,j,k,2), & - vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & + scon(i,j,k,3),tau(i,j,k),rho(i,j,k),lvel(i,j,k,1),lvel(i,j,k,2), & + lvel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det,x(i,j,k),y(i,j,k), & z(i,j,k),r(i,j,k),GRHydro_rho_min, GRHydro_reflevel, GRHydro_C2P_failed(i,j,k)) @@ -1600,18 +1600,18 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) @@ -1686,18 +1686,18 @@ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS) do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gxxl = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) @@ -1878,7 +1878,7 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS) write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) call CCTK_WARN(1,warnline) write(warnline,'(a20,3g16.7)') 'velocities: ',& - vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3) + lvel(i,j,k,1),lvel(i,j,k,2),lvel(i,j,k,3) call CCTK_WARN(0,warnline) !$OMP END CRITICAL diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90 index 0fcf71f..247b398 100644 --- a/src/GRHydro_ENOReconstruct_drv.F90 +++ b/src/GRHydro_ENOReconstruct_drv.F90 @@ -14,15 +14,15 @@ #include "SpaceMask.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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) +#define Bvecx(i,j,k) lBvec(i,j,k,1) +#define Bvecy(i,j,k) lBvec(i,j,k,2) +#define Bvecz(i,j,k) lBvec(i,j,k,3) #define Bconsx(i,j,k) Bcons(i,j,k,1) #define Bconsy(i,j,k) Bcons(i,j,k,2) #define Bconsz(i,j,k) Bcons(i,j,k,3) diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90 index 2c5776b..d1c022c 100644 --- a/src/GRHydro_EoSChangeGamma.F90 +++ b/src/GRHydro_EoSChangeGamma.F90 @@ -14,9 +14,9 @@ #include "cctk_Arguments.h" #include "GRHydro_Macros.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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) @@ -78,10 +78,10 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 prim2conpolytype(GRHydro_polytrope_handle,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 = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) + call prim2conpolytype(GRHydro_polytrope_handle,gaa(i,j,k),gab(i,j,k),& + gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(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)) @@ -242,10 +242,10 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 prim2con(GRHydro_eos_handle,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 = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) + call prim2con(GRHydro_eos_handle,gaa(i,j,k),gab(i,j,k),& + gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(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)) diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90 index 2074b60..ee6c0ad 100644 --- a/src/GRHydro_HLLE.F90 +++ b/src/GRHydro_HLLE.F90 @@ -104,26 +104,26 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) !!$ left and right points. if (flux_direction == 1) then - avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + & - betax(i,j,k)) + avg_beta = 0.5d0 * (betaa(i+xoffset,j+yoffset,k+zoffset) + & + betaa(i,j,k)) else if (flux_direction == 2) then - avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + & - betay(i,j,k)) + avg_beta = 0.5d0 * (betab(i+xoffset,j+yoffset,k+zoffset) + & + betab(i,j,k)) else if (flux_direction == 3) then - avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + & - betaz(i,j,k)) + avg_beta = 0.5d0 * (betac(i+xoffset,j+yoffset,k+zoffset) + & + betac(i,j,k)) else call CCTK_WARN(0, "Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) - gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k)) - gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k)) - gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k)) - gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k)) - gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) - gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) + gxxh = 0.5d0 * (gaa(i+xoffset,j+yoffset,k+zoffset) + gaa(i,j,k)) + gxyh = 0.5d0 * (gab(i+xoffset,j+yoffset,k+zoffset) + gab(i,j,k)) + gxzh = 0.5d0 * (gac(i+xoffset,j+yoffset,k+zoffset) + gac(i,j,k)) + gyyh = 0.5d0 * (gbb(i+xoffset,j+yoffset,k+zoffset) + gbb(i,j,k)) + gyzh = 0.5d0 * (gbc(i+xoffset,j+yoffset,k+zoffset) + gbc(i,j,k)) + gzzh = 0.5d0 * (gcc(i+xoffset,j+yoffset,k+zoffset) + gcc(i,j,k)) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\ gyyh,gyzh,gzzh) @@ -467,26 +467,26 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) !!$ left and right points. if (flux_direction == 1) then - avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + & - betax(i,j,k)) + avg_beta = 0.5d0 * (betaa(i+xoffset,j+yoffset,k+zoffset) + & + betaa(i,j,k)) else if (flux_direction == 2) then - avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + & - betay(i,j,k)) + avg_beta = 0.5d0 * (betab(i+xoffset,j+yoffset,k+zoffset) + & + betab(i,j,k)) else if (flux_direction == 3) then - avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + & - betaz(i,j,k)) + avg_beta = 0.5d0 * (betac(i+xoffset,j+yoffset,k+zoffset) + & + betac(i,j,k)) else call CCTK_WARN(0, "Flux direction not x,y,z") end if avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) - gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k)) - gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k)) - gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k)) - gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k)) - gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) - gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) + gxxh = 0.5d0 * (gaa(i+xoffset,j+yoffset,k+zoffset) + gaa(i,j,k)) + gxyh = 0.5d0 * (gab(i+xoffset,j+yoffset,k+zoffset) + gab(i,j,k)) + gxzh = 0.5d0 * (gac(i+xoffset,j+yoffset,k+zoffset) + gac(i,j,k)) + gyyh = 0.5d0 * (gbb(i+xoffset,j+yoffset,k+zoffset) + gbb(i,j,k)) + gyzh = 0.5d0 * (gbc(i+xoffset,j+yoffset,k+zoffset) + gbc(i,j,k)) + gzzh = 0.5d0 * (gcc(i+xoffset,j+yoffset,k+zoffset) + gcc(i,j,k)) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\ gyyh,gyzh,gzzh) diff --git a/src/GRHydro_Jacobian_state.c b/src/GRHydro_Jacobian_state.c new file mode 100644 index 0000000..4dd6f47 --- /dev/null +++ b/src/GRHydro_Jacobian_state.c @@ -0,0 +1,32 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + + + + +/** + Checks states for Jacobians +*/ +void +GRHydro_check_Jacobian_state (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + int ierr; + + if (*jacobian_state == 0) + { + CCTK_WARN(0, "GRHydro requires storage for Jacobians. Please tell your Coordinates implementation to store the Jacobians!"); + } + + if (*inverse_jacobian_state == 0) + { + CCTK_WARN(0, "GRHydro requires storage for inverse Jacobians. Please tell your Coordinates implementation to store the inverse Jacobians!"); + } + +} + + + diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90 index 963dbf6..a8058cc 100644 --- a/src/GRHydro_PPMMReconstruct_drv.F90 +++ b/src/GRHydro_PPMMReconstruct_drv.F90 @@ -14,15 +14,15 @@ #include "SpaceMask.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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) +#define Bvecx(i,j,k) lBvec(i,j,k,1) +#define Bvecy(i,j,k) lBvec(i,j,k,2) +#define Bvecz(i,j,k) lBvec(i,j,k,3) #define Bconsx(i,j,k) Bcons(i,j,k,1) #define Bconsy(i,j,k) Bcons(i,j,k,2) #define Bconsz(i,j,k) Bcons(i,j,k,3) @@ -153,7 +153,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),psidcplus(:,j,k),epsplus(:,j,k),& 1,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & - gzz(:,j,k), betax(:,j,k), alp(:,j,k),& + gzz(:,j,k), betaa(:,j,k), alp(:,j,k),& w_lorentz(:,j,k), & 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & GRHydro_mppm_eigenvalue_x_right, & @@ -168,7 +168,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) Bvecxplus(:,j,k),Bvecyplus(:,j,k),Bveczplus(:,j,k),dump(:,j,k),epsplus(:,j,k),& 0,trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & - gzz(:,j,k), betax(:,j,k), alp(:,j,k),& + gzz(:,j,k), betaa(:,j,k), alp(:,j,k),& w_lorentz(:,j,k), & 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & GRHydro_mppm_eigenvalue_x_right, & @@ -224,7 +224,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),psidcplus(j,:,k),epsplus(j,:,k),& 1,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & - gxx(j,:,k), betay(j,:,k), alp(j,:,k),& + gxx(j,:,k), betab(j,:,k), alp(j,:,k),& w_lorentz(j,:,k), & 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & GRHydro_mppm_eigenvalue_y_right, & @@ -239,7 +239,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) Bvecyplus(j,:,k),Bveczplus(j,:,k),Bvecxplus(j,:,k),dump(j,:,k),epsplus(j,:,k),& 0,trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & - gxx(j,:,k), betay(j,:,k), alp(j,:,k),& + gxx(j,:,k), betab(j,:,k), alp(j,:,k),& w_lorentz(j,:,k), & 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & GRHydro_mppm_eigenvalue_y_right, & @@ -294,7 +294,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),psidcplus(j,k,:),epsplus(j,k,:),& 1,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & - gyy(j,k,:), betaz(j,k,:), alp(j,k,:),& + gyy(j,k,:), betac(j,k,:), alp(j,k,:),& w_lorentz(j,k,:), & 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & GRHydro_mppm_eigenvalue_z_right, & @@ -309,7 +309,7 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) Bveczplus(j,k,:),Bvecxplus(j,k,:),Bvecyplus(j,k,:),dump(j,k,:),epsplus(j,k,:),& 0,trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & - gyy(j,k,:), betaz(j,k,:), alp(j,k,:),& + gyy(j,k,:), betac(j,k,:), alp(j,k,:),& w_lorentz(j,k,:), & 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & GRHydro_mppm_eigenvalue_z_right, & diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90 index b514247..2ecbb71 100644 --- a/src/GRHydro_PPMReconstruct_drv.F90 +++ b/src/GRHydro_PPMReconstruct_drv.F90 @@ -14,9 +14,9 @@ #include "SpaceMask.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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) @@ -125,8 +125,8 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) velzminus(:,j,k),epsminus(:,j,k),rhoplus(:,j,k),& velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),epsplus(:,j,k),& trivial_rp(:,j,k), hydro_excision_mask(:,j,k),& - gxx(:,j,k), gxy(:,j,k), gxz(:,j,k), gyy(:,j,k), gyz(:,j,k), & - gzz(:,j,k), betax(:,j,k), alp(:,j,k),& + gaa(:,j,k), gab(:,j,k), gac(:,j,k), gbb(:,j,k), gbc(:,j,k), & + gcc(:,j,k), betaa(:,j,k), alp(:,j,k),& w_lorentz(:,j,k), & 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & GRHydro_mppm_eigenvalue_x_right, & @@ -176,8 +176,8 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) velxminus(j,:,k),epsminus(j,:,k),rhoplus(j,:,k),& velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),epsplus(j,:,k),& trivial_rp(j,:,k), hydro_excision_mask(j,:,k),& - gyy(j,:,k), gyz(j,:,k), gxy(j,:,k), gzz(j,:,k), gxz(j,:,k), & - gxx(j,:,k), betay(j,:,k), alp(j,:,k),& + gbb(j,:,k), gbc(j,:,k), gab(j,:,k), gcc(j,:,k), gac(j,:,k), & + gaa(j,:,k), betab(j,:,k), alp(j,:,k),& w_lorentz(j,:,k), & 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & GRHydro_mppm_eigenvalue_y_right, & @@ -227,8 +227,8 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) velyminus(j,k,:),epsminus(j,k,:),rhoplus(j,k,:),& velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),epsplus(j,k,:),& trivial_rp(j,k,:), hydro_excision_mask(j,k,:),& - gzz(j,k,:), gxz(j,k,:), gyz(j,k,:), gxx(j,k,:), gxy(j,k,:), & - gyy(j,k,:), betaz(j,k,:), alp(j,k,:),& + gcc(j,k,:), gac(j,k,:), gbc(j,k,:), gaa(j,k,:), gab(j,k,:), & + gbb(j,k,:), betac(j,k,:), alp(j,k,:),& w_lorentz(j,k,:), & 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & GRHydro_mppm_eigenvalue_z_right, & diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index fbe118b..2bd2d29 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -13,9 +13,9 @@ #include "cctk_Functions.h" #include "GRHydro_Macros.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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) @@ -44,44 +44,44 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr CCTK_REAL :: xtemp(1) if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr,& - !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, & + !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) - call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,& - gyzl,gzzl, & + call prim2con(GRHydro_eos_handle, gaal,gabl,gacl,gbbl,& + gbcl,gccl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) - call prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + call prim2con(GRHydro_eos_handle, gaar,gabr,gacr,& + gbbr,gbcr,gccr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -94,27 +94,27 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) !$OMP END PARALLEL DO else !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& - !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, & + !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -124,8 +124,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i-xoffset,j-yoffset,k-zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxl,gxyl,gxzl,gyyl,& - gyzl,gzzl, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaal,gabl,gacl,gbbl,& + gbcl,gccl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & @@ -141,8 +141,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i+xoffset,j+yoffset,k+zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, & - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaar,gabr,gacr,& + gbbr,gbcr,gccr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -353,12 +353,12 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) - call prim2con(GRHydro_eos_handle,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 prim2con(GRHydro_eos_handle,gaa(i,j,k),& + gab(i,j,k),gac(i,j,k),& + gbb(i,j,k),gbc(i,j,k),gcc(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)) @@ -372,13 +372,13 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,i,j,k,& - x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),& - gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),& + gab(i,j,k),gac(i,j,k),& + gbb(i,j,k),gbc(i,j,k),gcc(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), & @@ -429,40 +429,40 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr - !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) - - call prim2conpolytype(GRHydro_eos_handle, gxxl,gxyl,gxzl,& - gyyl,gyzl,gzzl, & + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) + + call prim2conpolytype(GRHydro_eos_handle, gaal,gabl,gacl,& + gbbl,gbcl,gccl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) - call prim2conpolytype(GRHydro_eos_handle, gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + call prim2conpolytype(GRHydro_eos_handle, gaar,gabr,gacr,& + gbbr,gbcr,gccr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -586,11 +586,11 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) - call prim2conpolytype(GRHydro_eos_handle,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 prim2conpolytype(GRHydro_eos_handle,gaa(i,j,k),gab(i,j,k),& + gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(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)) @@ -637,46 +637,46 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * & sqrt(avg_detr) * rhoplus(i,j,k) / & sqrt(1.d0 - & - (gxxr * velxplus(i,j,k)**2 + & - gyyr * velyplus(i,j,k)**2 + & - gzzr * velzplus(i,j,k)**2 + & - 2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + & - gxzr * velxplus(i,j,k) * velzplus(i,j,k) + & - gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) + (gaar * velxplus(i,j,k)**2 + & + gbbr * velyplus(i,j,k)**2 + & + gccr * velzplus(i,j,k)**2 + & + 2.d0 * (gabr * velxplus(i,j,k) * velyplus(i,j,k) + & + gacr * velxplus(i,j,k) * velzplus(i,j,k) + & + gbcr * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * & sqrt(avg_detl) * rhominus(i,j,k) / & sqrt(1.d0 - & - (gxxl * velxminus(i,j,k)**2 + & - gyyl * velyminus(i,j,k)**2 + & - gzzl * velzminus(i,j,k)**2 + & - 2.d0 * (gxyl * velxminus(i,j,k) * velyminus(i,j,k) + & - gxzl * velxminus(i,j,k) * velzminus(i,j,k) + & - gyzl * velyminus(i,j,k) * velzminus(i,j,k) ) ) ) + (gaal * velxminus(i,j,k)**2 + & + gbbl * velyminus(i,j,k)**2 + & + gccl * velzminus(i,j,k)**2 + & + 2.d0 * (gabl * velxminus(i,j,k) * velyminus(i,j,k) + & + gacl * velxminus(i,j,k) * velzminus(i,j,k) + & + gbcl * velyminus(i,j,k) * velzminus(i,j,k) ) ) ) end do end do diff --git a/src/GRHydro_RegisterVars.cc b/src/GRHydro_RegisterVars.cc index 0d99955..5e91721 100644 --- a/src/GRHydro_RegisterVars.cc +++ b/src/GRHydro_RegisterVars.cc @@ -53,8 +53,13 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) register_constrained("HydroBase::rho"); register_constrained("HydroBase::press"); register_constrained("HydroBase::eps"); - register_constrained("HydroBase::vel"); register_constrained("HydroBase::w_lorentz"); + register_constrained("HydroBase::vel"); + register_constrained("GRHydro::lvel"); + + register_constrained("grhydro::local_shift"); + register_constrained("grhydro::local_metric"); + register_constrained("grhydro::local_extrinsic_curvature"); if (CCTK_EQUALS(evolution_method, "GRHydro")) { @@ -64,6 +69,7 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { register_constrained("HydroBase::Bvec"); + register_constrained("GRHydro::lBvec"); register_evolved("GRHydro::Bcons", "GRHydro::Bconsrhs"); if(clean_divergence) { register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs"); diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90 index f8fd9eb..969e88f 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -9,27 +9,27 @@ ! Second order f.d. -#define DIFF_X_2(q) 0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx -#define DIFF_Y_2(q) 0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy -#define DIFF_Z_2(q) 0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz +#define DIFF_X_2(q) 0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * ida +#define DIFF_Y_2(q) 0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idb +#define DIFF_Z_2(q) 0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idc ! Fourth order f.d. #define DIFF_X_4(q) (-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ - q(i-2,j,k)) / 12.d0 * idx + q(i-2,j,k)) / 12.d0 * ida #define DIFF_Y_4(q) (-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ - q(i,j-2,k)) / 12.d0 * idy + q(i,j-2,k)) / 12.d0 * idb #define DIFF_Z_4(q) (-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ - q(i,j,k-2)) / 12.d0 * idz + q(i,j,k-2)) / 12.d0 * idc #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "GRHydro_Macros.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 vela(i,j,k) lvel(i,j,k,1) +#define velb(i,j,k) lvel(i,j,k,2) +#define velc(i,j,k) lvel(i,j,k,3) /*@@ @routine SourceTerms @date Sat Jan 26 02:04:21 2002 @@ -52,26 +52,26 @@ subroutine SourceTerms(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - CCTK_INT :: i, j, k, nx, ny, nz + CCTK_INT :: i, j, k, na, nb, nc CCTK_REAL :: one, two, half - CCTK_REAL :: t00, t0x, t0y, t0z, txx, txy, txz, tyy, tyz, tzz - CCTK_REAL :: sqrtdet, det, uxx, uxy, uxz, uyy, uyz, uzz, rhoenthalpyW2 - CCTK_REAL :: shiftx, shifty, shiftz, velxshift, velyshift, velzshift - CCTK_REAL :: vlowx, vlowy, vlowz - CCTK_REAL :: dx_betax, dx_betay, dx_betaz, dy_betax, dy_betay,& - dy_betaz, dz_betax, dz_betay, dz_betaz - CCTK_REAL :: dx_alp, dy_alp, dz_alp - CCTK_REAL :: tau_source, sx_source, sy_source, sz_source - CCTK_REAL :: localgxx,localgxy,localgxz,localgyy,localgyz,localgzz - CCTK_REAL :: dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz - CCTK_REAL :: dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz - CCTK_REAL :: dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz - CCTK_REAL :: dx, dy, dz, idx, idy, idz - CCTK_REAL :: psi4pt, dx_psi4, dy_psi4, dz_psi4 - CCTK_REAL :: shiftshiftk, shiftkx, shiftky, shiftkz + CCTK_REAL :: t00, t0a, t0b, t0c, taa, tab, tac, tbb, tbc, tcc + CCTK_REAL :: sqrtdet, det, uaa, uab, uac, ubb, ubc, ucc, rhoenthalpyW2 + CCTK_REAL :: shifta, shiftb, shiftc, velashift, velbshift, velcshift + CCTK_REAL :: vlowa, vlowb, vlowc + CCTK_REAL :: da_betaa, da_betab, da_betac, db_betaa, db_betab,& + db_betac, dc_betaa, dc_betab, dc_betac + CCTK_REAL :: da_alp, db_alp, dc_alp + CCTK_REAL :: tau_source, sa_source, sb_source, sc_source + CCTK_REAL :: localgaa,localgab,localgac,localgbb,localgbc,localgcc + CCTK_REAL :: da_gaa, da_gab, da_gac, da_gbb, da_gbc, da_gcc + CCTK_REAL :: db_gaa, db_gab, db_gac, db_gbb, db_gbc, db_gcc + CCTK_REAL :: dc_gaa, dc_gab, dc_gac, dc_gbb, dc_gbc, dc_gcc + CCTK_REAL :: da, db, dc, ida, idb, idc + CCTK_REAL :: psi4pt, da_psi4, db_psi4, dc_psi4 + CCTK_REAL :: shiftshiftk, shiftka, shiftkb, shiftkc CCTK_REAL :: sumTK - CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz - CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz + CCTK_REAL :: halfshiftdga, halfshiftdgb, halfshiftdgc + CCTK_REAL :: halfTdga, halfTdgb, halfTdgc CCTK_REAL :: invalp, invalp2 logical, allocatable, dimension (:,:,:) :: force_spatial_second_order @@ -79,15 +79,15 @@ subroutine SourceTerms(CCTK_ARGUMENTS) one = 1.0d0 two = 2.0d0 half = 0.5d0 - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - dx = CCTK_DELTA_SPACE(1) - dy = CCTK_DELTA_SPACE(2) - dz = CCTK_DELTA_SPACE(3) - idx = 1.d0/dx - idy = 1.d0/dy - idz = 1.d0/dz + na = cctk_lsh(1) + nb = cctk_lsh(2) + nc = cctk_lsh(3) + da = CCTK_DELTA_SPACE(1) + db = CCTK_DELTA_SPACE(2) + dc = CCTK_DELTA_SPACE(3) + ida = 1.d0/da + idb = 1.d0/db + idc = 1.d0/dc !!$ Initialize the update terms to be zero. !!$ This will guarantee that no garbage in the boundaries is updated. @@ -108,14 +108,14 @@ subroutine SourceTerms(CCTK_ARGUMENTS) !!$ differencing at boundaries and near excision regions. !!$ Copied straight from BSSN. - allocate (force_spatial_second_order(nx,ny,nz)) + allocate (force_spatial_second_order(na,nb,nc)) force_spatial_second_order = .FALSE. if (spatial_order > 2) then !$OMP PARALLEL DO PRIVATE(i, j, k) - do k = 1 + GRHydro_stencil, nz - GRHydro_stencil - do j = 1 + GRHydro_stencil, ny - GRHydro_stencil - do i = 1 + GRHydro_stencil, nx - GRHydro_stencil + do k = 1 + GRHydro_stencil, nc - GRHydro_stencil + do j = 1 + GRHydro_stencil, nb - GRHydro_stencil + do i = 1 + GRHydro_stencil, na - GRHydro_stencil if ((i < 3).or.(i > cctk_lsh(1) - 2).or. & (j < 3).or.(j > cctk_lsh(2) - 2).or. & (k < 3).or.(k > cctk_lsh(3) - 2) ) then @@ -132,22 +132,22 @@ subroutine SourceTerms(CCTK_ARGUMENTS) end if !$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,& - !$OMP localgxx,localgxy,localgxz,localgyy,localgyz,localgzz,& - !$OMP psi4pt,det,sqrtdet,rhoenthalpyW2,shiftx,shifty,shiftz,& - !$OMP dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz,& - !$OMP dz_betax,dz_betay,dz_betaz,velxshift,velyshift,velzshift,& - !$OMP vlowx,vlowy,vlowz,t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz,& - !$OMP dx_alp,dy_alp,dz_alp,tau_source,sx_source,sy_source,sz_source,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz,& - !$OMP dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz,& - !$OMP dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz,& - !$OMP dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz,& - !$OMP dx_psi4,dy_psi4,dz_psi4,shiftshiftk,shiftkx,shiftky,shiftkz,& - !$OMP sumTK,halfshiftdgx,halfshiftdgy,halfshiftdgz,& - !$OMP halfTdgx,halfTdgy,halfTdgz,invalp,invalp2) - do k=1 + GRHydro_stencil,nz - GRHydro_stencil - do j=1 + GRHydro_stencil,ny - GRHydro_stencil - do i=1 + GRHydro_stencil,nx - GRHydro_stencil + !$OMP localgaa,localgab,localgac,localgbb,localgbc,localgcc,& + !$OMP psi4pt,det,sqrtdet,rhoenthalpyW2,shifta,shiftb,shiftc,& + !$OMP da_betaa,da_betab,da_betac,db_betaa,db_betab,db_betac,& + !$OMP dc_betaa,dc_betab,dc_betac,velashift,velbshift,velcshift,& + !$OMP vlowa,vlowb,vlowc,t00,t0a,t0b,t0c,taa,tab,tac,tbb,tbc,tcc,& + !$OMP da_alp,db_alp,dc_alp,tau_source,sa_source,sb_source,sc_source,& + !$OMP uaa, uab, uac, ubb, ubc, ucc,& + !$OMP da_gaa, da_gab, da_gac, da_gbb, da_gbc, da_gcc,& + !$OMP db_gaa, db_gab, db_gac, db_gbb, db_gbc, db_gcc,& + !$OMP dc_gaa, dc_gab, dc_gac, dc_gbb, dc_gbc, dc_gcc,& + !$OMP da_psi4,db_psi4,dc_psi4,shiftshiftk,shiftka,shiftkb,shiftkc,& + !$OMP sumTK,halfshiftdga,halfshiftdgb,halfshiftdgc,& + !$OMP halfTdga,halfTdgb,halfTdgc,invalp,invalp2) + do k=1 + GRHydro_stencil,nc - GRHydro_stencil + do j=1 + GRHydro_stencil,nb - GRHydro_stencil + do i=1 + GRHydro_stencil,na - GRHydro_stencil local_spatial_order = spatial_order if (force_spatial_second_order(i,j,k)) then @@ -156,18 +156,18 @@ subroutine SourceTerms(CCTK_ARGUMENTS) !!$ Set the metric terms. - localgxx = gxx(i,j,k) - localgxy = gxy(i,j,k) - localgxz = gxz(i,j,k) - localgyy = gyy(i,j,k) - localgyz = gyz(i,j,k) - localgzz = gzz(i,j,k) + localgaa = gaa(i,j,k) + localgab = gab(i,j,k) + localgac = gac(i,j,k) + localgbb = gbb(i,j,k) + localgbc = gbc(i,j,k) + localgcc = gcc(i,j,k) - det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\ - localgyy, localgyz, localgzz) + det = SPATIAL_DETERMINANT(localgaa, localgab, localgac,\ + localgbb, localgbc, localgcc) sqrtdet = sqrt(det) - call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,& - localgxy, localgxz, localgyy, localgyz, localgzz) + call UpperMetric(uaa, uab, uac, ubb, ubc, ucc, det, localgaa,& + localgab, localgac, localgbb, localgbc, localgcc) !!$ All the matter variables (except velocity) always appear !!$ together in this form @@ -175,218 +175,218 @@ subroutine SourceTerms(CCTK_ARGUMENTS) rhoenthalpyW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k))*& w_lorentz(i,j,k)**2 - shiftx = betax(i,j,k) - shifty = betay(i,j,k) - shiftz = betaz(i,j,k) + shifta = betaa(i,j,k) + shiftb = betab(i,j,k) + shiftc = betac(i,j,k) if (local_spatial_order .eq. 2) then - dx_betax = DIFF_X_2(betax) - dx_betay = DIFF_X_2(betay) - dx_betaz = DIFF_X_2(betaz) + da_betaa = DIFF_X_2(betaa) + da_betab = DIFF_X_2(betab) + da_betac = DIFF_X_2(betac) - dy_betax = DIFF_Y_2(betax) - dy_betay = DIFF_Y_2(betay) - dy_betaz = DIFF_Y_2(betaz) + db_betaa = DIFF_Y_2(betaa) + db_betab = DIFF_Y_2(betab) + db_betac = DIFF_Y_2(betac) - dz_betax = DIFF_Z_2(betax) - dz_betay = DIFF_Z_2(betay) - dz_betaz = DIFF_Z_2(betaz) + dc_betaa = DIFF_Z_2(betaa) + dc_betab = DIFF_Z_2(betab) + dc_betac = DIFF_Z_2(betac) else - dx_betax = DIFF_X_4(betax) - dx_betay = DIFF_X_4(betay) - dx_betaz = DIFF_X_4(betaz) + da_betaa = DIFF_X_4(betaa) + da_betab = DIFF_X_4(betab) + da_betac = DIFF_X_4(betac) - dy_betax = DIFF_Y_4(betax) - dy_betay = DIFF_Y_4(betay) - dy_betaz = DIFF_Y_4(betaz) + db_betaa = DIFF_Y_4(betaa) + db_betab = DIFF_Y_4(betab) + db_betac = DIFF_Y_4(betac) - dz_betax = DIFF_Z_4(betax) - dz_betay = DIFF_Z_4(betay) - dz_betaz = DIFF_Z_4(betaz) + dc_betaa = DIFF_Z_4(betaa) + dc_betab = DIFF_Z_4(betab) + dc_betac = DIFF_Z_4(betac) end if invalp = 1.0d0 / alp(i,j,k) invalp2 = invalp**2 - velxshift = velx(i,j,k) - shiftx*invalp - velyshift = vely(i,j,k) - shifty*invalp - velzshift = velz(i,j,k) - shiftz*invalp - vlowx = velx(i,j,k)*localgxx + vely(i,j,k)*localgxy +& - velz(i,j,k)*localgxz - vlowy = velx(i,j,k)*localgxy + vely(i,j,k)*localgyy +& - velz(i,j,k)*localgyz - vlowz = velx(i,j,k)*localgxz + vely(i,j,k)*localgyz +& - velz(i,j,k)*localgzz + velashift = vela(i,j,k) - shifta*invalp + velbshift = velb(i,j,k) - shiftb*invalp + velcshift = velc(i,j,k) - shiftc*invalp + vlowa = vela(i,j,k)*localgaa + velb(i,j,k)*localgab +& + velc(i,j,k)*localgac + vlowb = vela(i,j,k)*localgab + velb(i,j,k)*localgbb +& + velc(i,j,k)*localgbc + vlowc = vela(i,j,k)*localgac + velb(i,j,k)*localgbc +& + velc(i,j,k)*localgcc !!$ For a change, these are T^{ij} t00 = (rhoenthalpyW2 - press(i,j,k))*invalp2 - t0x = rhoenthalpyW2*velxshift/alp(i,j,k) +& - press(i,j,k)*shiftx*invalp2 - t0y = rhoenthalpyW2*velyshift/alp(i,j,k) +& - press(i,j,k)*shifty*invalp2 - t0z = rhoenthalpyW2*velzshift/alp(i,j,k) +& - press(i,j,k)*shiftz*invalp2 - txx = rhoenthalpyW2*velxshift*velxshift +& - press(i,j,k)*(uxx - shiftx*shiftx*invalp2) - txy = rhoenthalpyW2*velxshift*velyshift +& - press(i,j,k)*(uxy - shiftx*shifty*invalp2) - txz = rhoenthalpyW2*velxshift*velzshift +& - press(i,j,k)*(uxz - shiftx*shiftz*invalp2) - tyy = rhoenthalpyW2*velyshift*velyshift +& - press(i,j,k)*(uyy - shifty*shifty*invalp2) - tyz = rhoenthalpyW2*velyshift*velzshift +& - press(i,j,k)*(uyz - shifty*shiftz*invalp2) - tzz = rhoenthalpyW2*velzshift*velzshift +& - press(i,j,k)*(uzz - shiftz*shiftz*invalp2) + t0a = rhoenthalpyW2*velashift/alp(i,j,k) +& + press(i,j,k)*shifta*invalp2 + t0b = rhoenthalpyW2*velbshift/alp(i,j,k) +& + press(i,j,k)*shiftb*invalp2 + t0c = rhoenthalpyW2*velcshift/alp(i,j,k) +& + press(i,j,k)*shiftc*invalp2 + taa = rhoenthalpyW2*velashift*velashift +& + press(i,j,k)*(uaa - shifta*shifta*invalp2) + tab = rhoenthalpyW2*velashift*velbshift +& + press(i,j,k)*(uab - shifta*shiftb*invalp2) + tac = rhoenthalpyW2*velashift*velcshift +& + press(i,j,k)*(uac - shifta*shiftc*invalp2) + tbb = rhoenthalpyW2*velbshift*velbshift +& + press(i,j,k)*(ubb - shiftb*shiftb*invalp2) + tbc = rhoenthalpyW2*velbshift*velcshift +& + press(i,j,k)*(ubc - shiftb*shiftc*invalp2) + tcc = rhoenthalpyW2*velcshift*velcshift +& + press(i,j,k)*(ucc - shiftc*shiftc*invalp2) !!$ Derivatives of the lapse, metric and shift if (local_spatial_order .eq. 2) then - dx_alp = DIFF_X_2(alp) - dy_alp = DIFF_Y_2(alp) - dz_alp = DIFF_Z_2(alp) + da_alp = DIFF_X_2(alp) + db_alp = DIFF_Y_2(alp) + dc_alp = DIFF_Z_2(alp) else - dx_alp = DIFF_X_4(alp) - dy_alp = DIFF_Y_4(alp) - dz_alp = DIFF_Z_4(alp) + da_alp = DIFF_X_4(alp) + db_alp = DIFF_Y_4(alp) + dc_alp = DIFF_Z_4(alp) end if if (local_spatial_order .eq. 2) then - dx_gxx = DIFF_X_2(gxx) - dx_gxy = DIFF_X_2(gxy) - dx_gxz = DIFF_X_2(gxz) - dx_gyy = DIFF_X_2(gyy) - dx_gyz = DIFF_X_2(gyz) - dx_gzz = DIFF_X_2(gzz) - dy_gxx = DIFF_Y_2(gxx) - dy_gxy = DIFF_Y_2(gxy) - dy_gxz = DIFF_Y_2(gxz) - dy_gyy = DIFF_Y_2(gyy) - dy_gyz = DIFF_Y_2(gyz) - dy_gzz = DIFF_Y_2(gzz) - dz_gxx = DIFF_Z_2(gxx) - dz_gxy = DIFF_Z_2(gxy) - dz_gxz = DIFF_Z_2(gxz) - dz_gyy = DIFF_Z_2(gyy) - dz_gyz = DIFF_Z_2(gyz) - dz_gzz = DIFF_Z_2(gzz) + da_gaa = DIFF_X_2(gaa) + da_gab = DIFF_X_2(gab) + da_gac = DIFF_X_2(gac) + da_gbb = DIFF_X_2(gbb) + da_gbc = DIFF_X_2(gbc) + da_gcc = DIFF_X_2(gcc) + db_gaa = DIFF_Y_2(gaa) + db_gab = DIFF_Y_2(gab) + db_gac = DIFF_Y_2(gac) + db_gbb = DIFF_Y_2(gbb) + db_gbc = DIFF_Y_2(gbc) + db_gcc = DIFF_Y_2(gcc) + dc_gaa = DIFF_Z_2(gaa) + dc_gab = DIFF_Z_2(gab) + dc_gac = DIFF_Z_2(gac) + dc_gbb = DIFF_Z_2(gbb) + dc_gbc = DIFF_Z_2(gbc) + dc_gcc = DIFF_Z_2(gcc) else - dx_gxx = DIFF_X_4(gxx) - dx_gxy = DIFF_X_4(gxy) - dx_gxz = DIFF_X_4(gxz) - dx_gyy = DIFF_X_4(gyy) - dx_gyz = DIFF_X_4(gyz) - dx_gzz = DIFF_X_4(gzz) - dy_gxx = DIFF_Y_4(gxx) - dy_gxy = DIFF_Y_4(gxy) - dy_gxz = DIFF_Y_4(gxz) - dy_gyy = DIFF_Y_4(gyy) - dy_gyz = DIFF_Y_4(gyz) - dy_gzz = DIFF_Y_4(gzz) - dz_gxx = DIFF_Z_4(gxx) - dz_gxy = DIFF_Z_4(gxy) - dz_gxz = DIFF_Z_4(gxz) - dz_gyy = DIFF_Z_4(gyy) - dz_gyz = DIFF_Z_4(gyz) - dz_gzz = DIFF_Z_4(gzz) + da_gaa = DIFF_X_4(gaa) + da_gab = DIFF_X_4(gab) + da_gac = DIFF_X_4(gac) + da_gbb = DIFF_X_4(gbb) + da_gbc = DIFF_X_4(gbc) + da_gcc = DIFF_X_4(gcc) + db_gaa = DIFF_Y_4(gaa) + db_gab = DIFF_Y_4(gab) + db_gac = DIFF_Y_4(gac) + db_gbb = DIFF_Y_4(gbb) + db_gbc = DIFF_Y_4(gbc) + db_gcc = DIFF_Y_4(gcc) + dc_gaa = DIFF_Z_4(gaa) + dc_gab = DIFF_Z_4(gab) + dc_gac = DIFF_Z_4(gac) + dc_gbb = DIFF_Z_4(gbb) + dc_gbc = DIFF_Z_4(gbc) + dc_gcc = DIFF_Z_4(gcc) end if -!!$ Contract the shift with the extrinsic curvature +!!$ Contract the shift with the eatrinsic curvature - shiftshiftk = shiftx*shiftx*kxx(i,j,k) + & - shifty*shifty*kyy(i,j,k) + & - shiftz*shiftz*kzz(i,j,k) + & - two*(shiftx*shifty*kxy(i,j,k) + & - shiftx*shiftz*kxz(i,j,k) + & - shifty*shiftz*kyz(i,j,k)) + shiftshiftk = shifta*shifta*kaa(i,j,k) + & + shiftb*shiftb*kbb(i,j,k) + & + shiftc*shiftc*kcc(i,j,k) + & + two*(shifta*shiftb*kab(i,j,k) + & + shifta*shiftc*kac(i,j,k) + & + shiftb*shiftc*kbc(i,j,k)) - shiftkx = shiftx*kxx(i,j,k) + shifty*kxy(i,j,k) + shiftz*kxz(i,j,k) - shiftky = shiftx*kxy(i,j,k) + shifty*kyy(i,j,k) + shiftz*kyz(i,j,k) - shiftkz = shiftx*kxz(i,j,k) + shifty*kyz(i,j,k) + shiftz*kzz(i,j,k) + shiftka = shifta*kaa(i,j,k) + shiftb*kab(i,j,k) + shiftc*kac(i,j,k) + shiftkb = shifta*kab(i,j,k) + shiftb*kbb(i,j,k) + shiftc*kbc(i,j,k) + shiftkc = shifta*kac(i,j,k) + shiftb*kbc(i,j,k) + shiftc*kcc(i,j,k) !!$ Contract the matter terms with the extrinsic curvature - sumTK = txx*kxx(i,j,k) + tyy*kyy(i,j,k) + tzz*kzz(i,j,k) & - + two*(txy*kxy(i,j,k) + txz*kxz(i,j,k) + tyz*kyz(i,j,k)) + sumTK = taa*kaa(i,j,k) + tbb*kbb(i,j,k) + tcc*kcc(i,j,k) & + + two*(tab*kab(i,j,k) + tac*kac(i,j,k) + tbc*kbc(i,j,k)) !!$ Update term for tau tau_source = t00* & - (shiftshiftk - (shiftx*dx_alp + shifty*dy_alp + shiftz*dz_alp) )& - + t0x*(-dx_alp + two*shiftkx) & - + t0y*(-dy_alp + two*shiftky) & - + t0z*(-dz_alp + two*shiftkz) & + (shiftshiftk - (shifta*da_alp + shiftb*db_alp + shiftc*dc_alp) )& + + t0a*(-da_alp + two*shiftka) & + + t0b*(-db_alp + two*shiftkb) & + + t0c*(-dc_alp + two*shiftkc) & + sumTK -!!$ The following looks very little like the terms in the +!!$ The following looks verb little like the terms in the !!$ standard papers. Take a look in the ThornGuide to see why !!$ it is really the same thing. !!$ Contract the shift with derivatives of the metric - halfshiftdgx = half*(shiftx*shiftx*dx_gxx + & - shifty*shifty*dx_gyy + shiftz*shiftz*dx_gzz) + & - shiftx*shifty*dx_gxy + shiftx*shiftz*dx_gxz + & - shifty*shiftz*dx_gyz - halfshiftdgy = half*(shiftx*shiftx*dy_gxx + & - shifty*shifty*dy_gyy + shiftz*shiftz*dy_gzz) + & - shiftx*shifty*dy_gxy + shiftx*shiftz*dy_gxz + & - shifty*shiftz*dy_gyz - halfshiftdgz = half*(shiftx*shiftx*dz_gxx + & - shifty*shifty*dz_gyy + shiftz*shiftz*dz_gzz) + & - shiftx*shifty*dz_gxy + shiftx*shiftz*dz_gxz + & - shifty*shiftz*dz_gyz + halfshiftdga = half*(shifta*shifta*da_gaa + & + shiftb*shiftb*da_gbb + shiftc*shiftc*da_gcc) + & + shifta*shiftb*da_gab + shifta*shiftc*da_gac + & + shiftb*shiftc*da_gbc + halfshiftdgb = half*(shifta*shifta*db_gaa + & + shiftb*shiftb*db_gbb + shiftc*shiftc*db_gcc) + & + shifta*shiftb*db_gab + shifta*shiftc*db_gac + & + shiftb*shiftc*db_gbc + halfshiftdgc = half*(shifta*shifta*dc_gaa + & + shiftb*shiftb*dc_gbb + shiftc*shiftc*dc_gcc) + & + shifta*shiftb*dc_gab + shifta*shiftc*dc_gac + & + shiftb*shiftc*dc_gbc !!$ Contract the matter with derivatives of the metric - halfTdgx = half*(txx*dx_gxx + tyy*dx_gyy + tzz*dx_gzz) +& - txy*dx_gxy + txz*dx_gxz + tyz*dx_gyz - halfTdgy = half*(txx*dy_gxx + tyy*dy_gyy + tzz*dy_gzz) +& - txy*dy_gxy + txz*dy_gxz + tyz*dy_gyz - halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) +& - txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz - - sx_source = t00*& - (halfshiftdgx - alp(i,j,k)*dx_alp) +& - t0x*(shiftx*dx_gxx + shifty*dx_gxy + shiftz*dx_gxz) +& - t0y*(shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) +& - t0z*(shiftx*dx_gxz + shifty*dx_gyz + shiftz*dx_gzz) +& - halfTdgx + rhoenthalpyW2*& - (vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz)*& + halfTdga = half*(taa*da_gaa + tbb*da_gbb + tcc*da_gcc) +& + tab*da_gab + tac*da_gac + tbc*da_gbc + halfTdgb = half*(taa*db_gaa + tbb*db_gbb + tcc*db_gcc) +& + tab*db_gab + tac*db_gac + tbc*db_gbc + halfTdgc = half*(taa*dc_gaa + tbb*dc_gbb + tcc*dc_gcc) +& + tab*dc_gab + tac*dc_gac + tbc*dc_gbc + + sa_source = t00*& + (halfshiftdga - alp(i,j,k)*da_alp) +& + t0a*(shifta*da_gaa + shiftb*da_gab + shiftc*da_gac) +& + t0b*(shifta*da_gab + shiftb*da_gbb + shiftc*da_gbc) +& + t0c*(shifta*da_gac + shiftb*da_gbc + shiftc*da_gcc) +& + halfTdga + rhoenthalpyW2*& + (vlowa*da_betaa + vlowb*da_betab + vlowc*da_betac)*& invalp - sy_source = t00*& - (halfshiftdgy - alp(i,j,k)*dy_alp) +& - t0x*(shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) +& - t0y*(shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) +& - t0z*(shiftx*dy_gxz + shifty*dy_gyz + shiftz*dy_gzz) +& - halfTdgy + rhoenthalpyW2*& - (vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz)*& + sb_source = t00*& + (halfshiftdgb - alp(i,j,k)*db_alp) +& + t0a*(shifta*db_gaa + shiftb*db_gab + shiftc*db_gac) +& + t0b*(shifta*db_gab + shiftb*db_gbb + shiftc*db_gbc) +& + t0c*(shifta*db_gac + shiftb*db_gbc + shiftc*db_gcc) +& + halfTdgb + rhoenthalpyW2*& + (vlowa*db_betaa + vlowb*db_betab + vlowc*db_betac)*& invalp - sz_source = t00*& - (halfshiftdgz - alp(i,j,k)*dz_alp) +& - t0x*(shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) +& - t0y*(shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) +& - t0z*(shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) +& - halfTdgz + rhoenthalpyW2*& - (vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz)*& + sc_source = t00*& + (halfshiftdgc - alp(i,j,k)*dc_alp) +& + t0a*(shifta*dc_gaa + shiftb*dc_gab + shiftc*dc_gac) +& + t0b*(shifta*dc_gab + shiftb*dc_gbb + shiftc*dc_gbc) +& + t0c*(shifta*dc_gac + shiftb*dc_gbc + shiftc*dc_gcc) +& + halfTdgc + rhoenthalpyW2*& + (vlowa*dc_betaa + vlowb*dc_betab + vlowc*dc_betac)*& invalp densrhs(i,j,k) = 0.d0 - srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source - srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source - srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source + srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sa_source + srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sb_source + srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sc_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source enddo diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90 index b314a62..96a71db 100644 --- a/src/GRHydro_TVDReconstruct_drv.F90 +++ b/src/GRHydro_TVDReconstruct_drv.F90 @@ -14,15 +14,15 @@ #include "SpaceMask.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 velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(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) +#define Bvecx(i,j,k) lBvec(i,j,k,1) +#define Bvecy(i,j,k) lBvec(i,j,k,2) +#define Bvecz(i,j,k) lBvec(i,j,k,3) #define Bconsx(i,j,k) Bcons(i,j,k,1) #define Bconsy(i,j,k) Bcons(i,j,k,2) #define Bconsz(i,j,k) Bcons(i,j,k,3) @@ -155,20 +155,20 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + lvel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + lvel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - vel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + lvel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & eps, epsplus, epsminus, trivial_rp, hydro_excision_mask) if(evolve_mhd.ne.0) then call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) + lBvec(:,:,:,1), Bvecxplus, Bvecxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask) + lBvec(:,:,:,2), Bvecyplus, Bvecyminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - Bvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) + lBvec(:,:,:,3), Bveczplus, Bveczminus, trivial_rp, hydro_excision_mask) endif else if (CCTK_EQUALS(recon_vars,"conservative")) then diff --git a/src/GRHydro_TransformTensorBasis.F90 b/src/GRHydro_TransformTensorBasis.F90 new file mode 100644 index 0000000..1aeda46 --- /dev/null +++ b/src/GRHydro_TransformTensorBasis.F90 @@ -0,0 +1,211 @@ +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + + +subroutine GRHydroTransformADMToLocalBasis(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i,j,k + + if (jacobian_state .eq. 0) then + call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!") + end if + + if (inverse_jacobian_state .eq. 0) then + call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!") + end if + + !$OMP PARALLEL DO PRIVATE(i,j,k) + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + gaa(i,j,k) = iJ11(i,j,k)**2 * gxx(i,j,k) + & + 2.0d0 * iJ11(i,j,k) * iJ21(i,j,k) * gxy(i,j,k) + & + 2.0d0 * iJ11(i,j,k) * iJ31(i,j,k) * gxz(i,j,k) + iJ21(i,j,k)**2 * gyy(i,j,k) + & + 2.0d0 * iJ21(i,j,k) * iJ31(i,j,k) * gyz(i,j,k) + iJ31(i,j,k)**2 * gzz(i,j,k) + + gab(i,j,k) = iJ11(i,j,k) * iJ12(i,j,k) * gxx(i,j,k) + & + (iJ11(i,j,k) * iJ22(i,j,k) + iJ21(i,j,k) * iJ12(i,j,k)) * gxy(i,j,k) + & + (iJ11(i,j,k) * iJ32(i,j,k) + iJ31(i,j,k) * iJ12(i,j,k)) * gxz(i,j,k) + & + iJ21(i,j,k) * iJ22(i,j,k) * gyy(i,j,k) + (iJ21(i,j,k) * iJ32(i,j,k) + & + iJ31(i,j,k) * iJ22(i,j,k)) * gyz(i,j,k) + iJ31(i,j,k) * iJ32(i,j,k) * gzz(i,j,k) + + gac(i,j,k) = iJ11(i,j,k) * iJ13(i,j,k) * gxx(i,j,k) + (iJ11(i,j,k) * iJ23(i,j,k) + & + iJ21(i,j,k) * iJ13(i,j,k)) * gxy(i,j,k) + & + (iJ11(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ13(i,j,k)) * gxz(i,j,k) + & + iJ21(i,j,k) * iJ23(i,j,k) * gyy(i,j,k) + & + (iJ21(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ23(i,j,k)) * gyz(i,j,k) + & + iJ31(i,j,k) * iJ33(i,j,k) * gzz(i,j,k) + + gbb(i,j,k) = iJ12(i,j,k)**2 * gxx(i,j,k) + & + 2.0d0 * iJ12(i,j,k) * iJ22(i,j,k) * gxy(i,j,k) + & + 2.0d0 * iJ12(i,j,k) * iJ32(i,j,k) * gxz(i,j,k) + iJ22(i,j,k)**2 * gyy(i,j,k) + & + 2.0d0 * iJ22(i,j,k) * iJ32(i,j,k) * gyz(i,j,k) + iJ32(i,j,k)**2 * gzz(i,j,k) + + gbc(i,j,k) = iJ12(i,j,k) * iJ13(i,j,k) * gxx(i,j,k) + & + (iJ12(i,j,k) * iJ23(i,j,k) + iJ22(i,j,k) * iJ13(i,j,k)) * gxy(i,j,k) + & + (iJ12(i,j,k) * iJ33(i,j,k) + iJ32(i,j,k) * iJ13(i,j,k)) * gxz(i,j,k) + & + iJ22(i,j,k) * iJ23(i,j,k) * gyy(i,j,k) + & + (iJ22(i,j,k)* iJ33(i,j,k) + iJ32(i,j,k) * iJ23(i,j,k)) * gyz(i,j,k) + & + iJ32(i,j,k) * iJ33(i,j,k) * gzz(i,j,k) + + gcc(i,j,k) = iJ13(i,j,k)**2 * gxx(i,j,k) + & + 2.0d0 * iJ13(i,j,k) * iJ23(i,j,k) * gxy(i,j,k) + & + 2.0d0 * iJ13(i,j,k) * iJ33(i,j,k) * gxz(i,j,k) + & + iJ23(i,j,k)**2 * gyy(i,j,k) + & + 2.0d0 * iJ23(i,j,k) * iJ33(i,j,k) * gyz(i,j,k) + & + iJ33(i,j,k)**2 * gzz(i,j,k) + + ! Transform extrinsic curvature from global to local basis. + ! Since extrinsic has covariant indices, use inverse Jacobian + kaa(i,j,k) = iJ11(i,j,k)**2 * kxx(i,j,k) + & + 2.0d0 * iJ11(i,j,k) * iJ21(i,j,k) * kxy(i,j,k) + & + 2.0d0 * iJ11(i,j,k) * iJ31(i,j,k) * kxz(i,j,k) + & + iJ21(i,j,k)**2 * kyy(i,j,k) + & + 2.0d0 * iJ21(i,j,k) * iJ31(i,j,k) * kyz(i,j,k) + & + iJ31(i,j,k)**2 * kzz(i,j,k) + + kab(i,j,k) = iJ11(i,j,k) * iJ12(i,j,k) * kxx(i,j,k) + & + (iJ11(i,j,k) * iJ22(i,j,k) + iJ21(i,j,k) * iJ12(i,j,k)) * kxy(i,j,k) + & + (iJ11(i,j,k) * iJ32(i,j,k) + iJ31(i,j,k) * iJ12(i,j,k)) * kxz(i,j,k) + & + iJ21(i,j,k) * iJ22(i,j,k) * kyy(i,j,k) + & + (iJ21(i,j,k) * iJ32(i,j,k) + iJ31(i,j,k) * iJ22(i,j,k)) * kyz(i,j,k) + & + iJ31(i,j,k) * iJ32(i,j,k) * kzz(i,j,k) + + + kac(i,j,k) = iJ11(i,j,k) * iJ13(i,j,k) * kxx(i,j,k) + & + (iJ11(i,j,k) * iJ23(i,j,k) + iJ21(i,j,k) * iJ13(i,j,k)) * kxy(i,j,k) + & + (iJ11(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ13(i,j,k)) * kxz(i,j,k) + & + iJ21(i,j,k) * iJ23(i,j,k) * kyy(i,j,k) + & + (iJ21(i,j,k) * iJ33(i,j,k) + iJ31(i,j,k) * iJ23(i,j,k)) * kyz(i,j,k) + & + iJ31(i,j,k) * iJ33(i,j,k) * kzz(i,j,k) + + kbb(i,j,k) = iJ12(i,j,k)**2 * kxx(i,j,k) + & + 2.0d0 * iJ12(i,j,k) * iJ22(i,j,k) * kxy(i,j,k) + & + 2.0d0 * iJ12(i,j,k) * iJ32(i,j,k) * kxz(i,j,k) + & + iJ22(i,j,k)**2 * kyy(i,j,k) + & + 2.0d0 * iJ22(i,j,k) * iJ32(i,j,k) * kyz(i,j,k) + & + iJ32(i,j,k)**2 * kzz(i,j,k) + + kbc(i,j,k) = iJ12(i,j,k) * iJ13(i,j,k) * kxx(i,j,k) + & + (iJ12(i,j,k) * iJ23(i,j,k) + iJ22(i,j,k) * iJ13(i,j,k)) * kxy(i,j,k) + & + (iJ12(i,j,k) * iJ33(i,j,k) + iJ32(i,j,k) * iJ13(i,j,k)) * kxz(i,j,k) + & + iJ22(i,j,k) * iJ23(i,j,k) * kyy(i,j,k) + & + (iJ22(i,j,k) * iJ33(i,j,k) + iJ32(i,j,k) * iJ23(i,j,k)) * kyz(i,j,k) + & + iJ32(i,j,k) * iJ33(i,j,k) * kzz(i,j,k) + + kcc(i,j,k) = iJ13(i,j,k)**2 * kxx(i,j,k) + & + 2.0d0 * iJ13(i,j,k) * iJ23(i,j,k) * kxy(i,j,k) + & + 2.0d0 * iJ13(i,j,k) * iJ33(i,j,k) * kxz(i,j,k) + & + iJ23(i,j,k)**2 * kyy(i,j,k) + & + 2.0d0 * iJ23(i,j,k) * iJ33(i,j,k) * kyz(i,j,k) + & + iJ33(i,j,k)**2 * kzz(i,j,k) + + ! Transform shift from global to local basis. + ! Since shift has contravariant index, use Jacobian + betaa(i,j,k) = betax(i,j,k)*J11(i,j,k) + & + betay(i,j,k)*J12(i,j,k) + & + betaz(i,j,k)*J13(i,j,k) + betab(i,j,k) = betax(i,j,k)*J21(i,j,k) + & + betay(i,j,k)*J22(i,j,k) + & + betaz(i,j,k)*J23(i,j,k) + betac(i,j,k) = betax(i,j,k)*J31(i,j,k) + & + betay(i,j,k)*J32(i,j,k) + & + betaz(i,j,k)*J33(i,j,k) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end subroutine GRHydroTransformADMToLocalBasis + + + +subroutine GRHydroTransformPrimToLocalBasis(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i,j,k + + if (jacobian_state .eq. 0) then + call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!") + end if + + if (inverse_jacobian_state .eq. 0) then + call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!") + end if + + !call CCTK_INFO("Transforming primitives to local basis.") + + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + ! Transform primitive velocity from global to local basis. + ! Since velocity has contravariant index, use Jacobian + lvel(i,j,k,1) = vel(i,j,k,1)*J11(i,j,k) + vel(i,j,k,2)*J12(i,j,k) + vel(i,j,k,3)*J13(i,j,k) + lvel(i,j,k,2) = vel(i,j,k,1)*J21(i,j,k) + vel(i,j,k,2)*J22(i,j,k) + vel(i,j,k,3)*J23(i,j,k) + lvel(i,j,k,3) = vel(i,j,k,1)*J31(i,j,k) + vel(i,j,k,2)*J32(i,j,k) + vel(i,j,k,3)*J33(i,j,k) + + ! Transform primitive B-field from global to local basis. + ! Since B-field has contravariant index, use Jacobian +! lBvec(i,j,k,1) = Bvec(i,j,k,1)*J11(i,j,k) + Bvec(i,j,k,2)*J12(i,j,k) + Bvec(i,j,k,3)*J13(i,j,k) +! lBvec(i,j,k,2) = Bvec(i,j,k,1)*J21(i,j,k) + Bvec(i,j,k,2)*J22(i,j,k) + Bvec(i,j,k,3)*J23(i,j,k) +! lBvec(i,j,k,3) = Bvec(i,j,k,1)*J31(i,j,k) + Bvec(i,j,k,2)*J32(i,j,k) + Bvec(i,j,k,3)*J33(i,j,k) + end do + end do + end do + !$OMP END PARALLEL DO + +end subroutine GRHydroTransformPrimToLocalBasis + + + +subroutine GRHydroTransformPrimToGlobalBasis(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i,j,k + + if (jacobian_state .eq. 0) then + call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!") + end if + + if (inverse_jacobian_state .eq. 0) then + call CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!") + end if + + !call CCTK_INFO("Transforming primitives to global basis.") + + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + ! Transform primitive velocity from local to global basis. + ! Since velocity has contravariant index, use inverse Jacobian + vel(i,j,k,1) = lvel(i,j,k,1)*iJ11(i,j,k) + lvel(i,j,k,2)*iJ12(i,j,k) + lvel(i,j,k,3)*iJ13(i,j,k) + vel(i,j,k,2) = lvel(i,j,k,1)*iJ21(i,j,k) + lvel(i,j,k,2)*iJ22(i,j,k) + lvel(i,j,k,3)*iJ23(i,j,k) + vel(i,j,k,3) = lvel(i,j,k,1)*iJ31(i,j,k) + lvel(i,j,k,2)*iJ32(i,j,k) + lvel(i,j,k,3)*iJ33(i,j,k) + + end do + end do + end do + !$OMP END PARALLEL DO + +end subroutine GRHydroTransformPrimToGlobalBasis + + + + + diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90 index d6bf540..0bb092c 100644 --- a/src/GRHydro_UpdateMask.F90 +++ b/src/GRHydro_UpdateMask.F90 @@ -14,15 +14,15 @@ #include "GRHydro_Macros.h" #include "SpaceMask.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 velx_p(i,j,k) vel_p(i,j,k,1) -#define vely_p(i,j,k) vel_p(i,j,k,2) -#define velz_p(i,j,k) vel_p(i,j,k,3) -#define velx_p_p(i,j,k) vel_p_p(i,j,k,1) -#define vely_p_p(i,j,k) vel_p_p(i,j,k,2) -#define velz_p_p(i,j,k) vel_p_p(i,j,k,3) +#define velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(i,j,k,3) +#define velx_p(i,j,k) lvel_p(i,j,k,1) +#define vely_p(i,j,k) lvel_p(i,j,k,2) +#define velz_p(i,j,k) lvel_p(i,j,k,3) +#define velx_p_p(i,j,k) lvel_p_p(i,j,k,1) +#define vely_p_p(i,j,k) lvel_p_p(i,j,k,2) +#define velz_p_p(i,j,k) lvel_p_p(i,j,k,3) subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS) @@ -246,9 +246,9 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.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)) - + det = SPATIAL_DETERMINANT(gaa(i,j,k), gab(i,j,k), gac(i,j,k), \ + gbb(i,j,k), gbc(i,j,k), gcc(i,j,k)) + if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low temperature(i,j,k) = grhydro_hot_atmo_temp @@ -259,8 +259,8 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) press(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),gxy(i,j,k),& - gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),gab(i,j,k),& + gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k), & det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & 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),& @@ -268,19 +268,19 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k) else call prim2conpolytype(GRHydro_polytrope_handle, & - 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, & + gaa(i,j,k), gab(i,j,k), gac(i,j,k), & + gbb(i,j,k), gbc(i,j,k), gcc(i,j,k), det, & dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & 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)) if (wk_atmosphere .eq. 0) then - atmosphere_mask(i, j, k) = 0 - SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\ - not_atmosphere) + atmosphere_mask(i, j, k) = 0 + SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\ + not_atmosphere) end if endif - end if + end if end do end do @@ -330,8 +330,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) if (rho(i,j,k) .le. rho_min) then rho(i,j,k) = rho_min - 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)) + det = SPATIAL_DETERMINANT(gaa(i,j,k), gab(i,j,k), gac(i,j,k), \ + gbb(i,j,k), gbc(i,j,k), gcc(i,j,k)) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 @@ -347,8 +347,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) press(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),gxy(i,j,k),& - gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),gab(i,j,k),& + gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k), & det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & 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),& @@ -360,8 +360,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) call prim2conpolytype(eos_handle, & - 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, & + gaa(i,j,k), gab(i,j,k), gac(i,j,k), & + gbb(i,j,k), gbc(i,j,k), gcc(i,j,k), det, & dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & 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)) @@ -374,8 +374,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0 - det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \ - gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k)) + det = SPATIAL_DETERMINANT(gaa_p(i,j,k), gab_p(i,j,k), gac_p(i,j,k), \ + gbb_p(i,j,k), gbc_p(i,j,k), gcc_p(i,j,k)) if(evolve_temper.ne.0) then ! set the temperature to be relatively low @@ -387,8 +387,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) press_p(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p(i,j,k),gxy_p(i,j,k),& - gxz_p(i,j,k),gyy_p(i,j,k),gyz_p(i,j,k),gzz_p(i,j,k), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa_p(i,j,k),gab_p(i,j,k),& + gac_p(i,j,k),gbb_p(i,j,k),gbc_p(i,j,k),gcc_p(i,j,k), & det,dens_p(i,j,k),scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k),& @@ -400,8 +400,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) call prim2conpolytype(eos_handle, & - gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), & - gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, & + gaa_p(i,j,k), gab_p(i,j,k), gac_p(i,j,k), & + gbb_p(i,j,k), gbc_p(i,j,k), gcc_p(i,j,k), det, & dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k)) @@ -415,8 +415,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) velx_p_p(i,j,k) = 0.0d0 vely_p_p(i,j,k) = 0.0d0 velz_p_p(i,j,k) = 0.0d0 - det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \ - gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k)) + det = SPATIAL_DETERMINANT(gaa_p_p(i,j,k), gab_p_p(i,j,k), gac_p_p(i,j,k), \ + gbb_p_p(i,j,k), gbc_p_p(i,j,k), gcc_p_p(i,j,k)) if(evolve_temper.ne.0) then ! set the temperature to be relatively low @@ -427,8 +427,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),& press_p_p(i,j,k),keyerr,anyerr) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p_p(i,j,k),gxy_p_p(i,j,k),& - gxz_p_p(i,j,k),gyy_p_p(i,j,k),gyz_p_p(i,j,k),gzz_p_p(i,j,k), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaa_p_p(i,j,k),gab_p_p(i,j,k),& + gac_p_p(i,j,k),gbb_p_p(i,j,k),gbc_p_p(i,j,k),gcc_p_p(i,j,k), & det,dens_p_p(i,j,k),scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k),& @@ -440,8 +440,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) call prim2conpolytype(eos_handle, & - gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), & - gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, & + gaa_p_p(i,j,k), gab_p_p(i,j,k), gac_p_p(i,j,k), & + gbb_p_p(i,j,k), gbc_p_p(i,j,k), gcc_p_p(i,j,k), det, & dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k)) diff --git a/src/make.code.defn b/src/make.code.defn index a4e0873..6ff530b 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -10,7 +10,6 @@ SRCS = Utils.F90 \ GRHydro_DivergenceClean.F90 \ GRHydro_Eigenproblem.F90 \ GRHydro_Eigenproblem_Marquina.F90 \ - GRHydro_ENOReconstruct_drv.F90 \ GRHydro_ENOReconstruct.F90 \ GRHydro_ENOScalars.F90 \ GRHydro_EOS.c \ @@ -20,7 +19,6 @@ SRCS = Utils.F90 \ GRHydro_Loop.F90 \ GRHydro_Minima.F90 \ GRHydro_Minima.cc \ - GRHydro_PPMReconstruct_drv.F90 \ GRHydro_PPM.F90 \ GRHydro_ParamCheck.F90 \ GRHydro_Particle.F90 \ @@ -38,7 +36,6 @@ SRCS = Utils.F90 \ GRHydro_SlopeLimiter.F90 \ GRHydro_Source.F90 \ GRHydro_Startup.F90 \ - GRHydro_TVDReconstruct_drv.F90 \ GRHydro_TVDReconstruct.F90 \ GRHydro_Marquina.F90 \ GRHydro_UpdateMask.F90 \ @@ -52,15 +49,19 @@ SRCS = Utils.F90 \ GRHydro_EigenproblemM.F90 \ GRHydro_FluxM.F90 \ GRHydro_HLLEM.F90 \ - GRHydro_PPMMReconstruct_drv.F90 \ GRHydro_PPMM.F90 \ GRHydro_Prim2ConM.F90 \ GRHydro_RiemannSolveM.F90 \ GRHydro_SourceM.F90 \ GRHydro_TmunuM.F90 \ GRHydro_UpdateMaskM.F90 \ - GRHydro_UtilsM.F90 - + GRHydro_UtilsM.F90 \ + GRHydro_TransformTensorBasis.F90 \ + GRHydro_Jacobian_state.c \ + GRHydro_PPMMReconstruct_drv.F90\ + GRHydro_ENOReconstruct_drv.F90\ + GRHydro_PPMReconstruct_drv.F90\ + GRHydro_TVDReconstruct_drv.F90 -- cgit v1.2.3