diff options
| -rw-r--r-- | interface.ccl | 2 | ||||
| -rw-r--r-- | schedule.ccl | 81 | ||||
| -rw-r--r-- | src/GRHydro_Boundaries.F90 | 51 | ||||
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 244 | ||||
| -rw-r--r-- | src/GRHydro_EoSChangeGamma.F90 | 79 | ||||
| -rw-r--r-- | src/GRHydro_HLLE.F90 | 110 | ||||
| -rw-r--r-- | src/GRHydro_Jacobian_state.c | 22 | ||||
| -rw-r--r-- | src/GRHydro_PPMMReconstruct_drv.F90 | 65 | ||||
| -rw-r--r-- | src/GRHydro_PPMReconstruct_drv.F90 | 53 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheck.F90 | 9 | ||||
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 350 | ||||
| -rw-r--r-- | src/GRHydro_RegisterVars.cc | 18 | ||||
| -rw-r--r-- | src/GRHydro_Source.F90 | 167 | ||||
| -rw-r--r-- | src/GRHydro_TVDReconstruct_drv.F90 | 23 | ||||
| -rw-r--r-- | src/GRHydro_TransformTensorBasis.F90 | 211 | ||||
| -rw-r--r-- | src/GRHydro_TransformTensorBasis.c | 438 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMask.F90 | 175 | ||||
| -rw-r--r-- | src/make.code.defn | 2 |
18 files changed, 1503 insertions, 597 deletions
diff --git a/interface.ccl b/interface.ccl index 3a347bd..1fc2731 100644 --- a/interface.ccl +++ b/interface.ccl @@ -7,7 +7,7 @@ #################################################################### implements: GRHydro -inherits: ADMBase, Boundary, SpaceMask, ADMMacros, Tmunubase, HydroBase, Coordinates +inherits: ADMBase, Boundary, SpaceMask, ADMMacros, Tmunubase, HydroBase #INCLUDES SOURCE: GRHydro_EMTensor.inc in CalcTmunu.inc #INCLUDES: GRHydro_EMTensor_temps.inc in CalcTmunu_temps.inc diff --git a/schedule.ccl b/schedule.ccl index 1c64d4f..208648b 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -184,13 +184,20 @@ if (!CCTK_Equals(initial_shift,"none")) ### Storage for local tensor quantities ### ############################################## -STORAGE: lvel[3] -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { - STORAGE: lBvec[3] +# the official test is to test Coordinates::general_coordinates, however +# at schedule time, this is not yet set or accessible if we really insist on +# running with Coordinates but with a Cartesian (trivial) coordinate system +# _and_ we really want to safe this memory then it has to be turned off at +# runtime via CCTK_GroupStorageDecrease +if(CCTK_IsImplementationActive("Coordinates")) { + STORAGE: lvel[3] + if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { + STORAGE: lBvec[3] + } + STORAGE: local_metric[3] + STORAGE: local_extrinsic_curvature + STORAGE: local_shift } -STORAGE: local_metric[3] -STORAGE: local_extrinsic_curvature -STORAGE: local_shift ############################################## ### Storage for the conformal state scalar ### @@ -226,11 +233,13 @@ schedule GRHydro_ParamCheck AT PARAMCHECK LANG: Fortran } "Check parameters" -schedule GRHydro_check_Jacobian_state AT BASEGRID AFTER (TmunuBase_SetStressEnergyState Coordinates_SetGlobalCoords_Group) -{ - LANG: C - OPTIONS: GLOBAL -} "Test state of Jacobians" +if(CCTK_IsImplementationActive("Coordinates")) { + schedule GRHydro_check_Jacobian_state AT BASEGRID AFTER (TmunuBase_SetStressEnergyState Coordinates_SetGlobalCoords_Group) + { + LANG: C + OPTIONS: GLOBAL + } "Test state of Jacobians" +} ###################################### ### Standard symmetry registration ### @@ -501,39 +510,41 @@ schedule GRHydro_SetupDescriptors AT CCTK_Initial BEFORE HydroBase_Initial ##################################################################### -schedule GRHydroTransformPrimToLocalBasis AT INITIAL AFTER (HydroBase_Initial, ADMBase_PostInitial) BEFORE HydroBase_Prim2ConInitial -{ - LANG: FORTRAN -} "Transform primitive vars to local tensor basis." +if(CCTK_IsImplementationActive("Coordinates")) { + schedule GRHydroTransformPrimToLocalBasis AT INITIAL AFTER (HydroBase_Initial, ADMBase_PostInitial) BEFORE HydroBase_Prim2ConInitial + { + LANG: C + } "Transform primitive vars to local tensor basis." -schedule GRHydroTransformADMToLocalBasis AT INITIAL AFTER HydroBase_Initial BEFORE GRHydroTransformPrimToLocalBasis -{ - LANG: FORTRAN -} "Transform ADM metric, extr. curv. and shift to local tensor basis." + schedule GRHydroTransformADMToLocalBasis AT INITIAL AFTER HydroBase_Initial BEFORE GRHydroTransformPrimToLocalBasis + { + LANG: C + } "Transform ADM metric, extr. curv. and shift to local tensor basis." -schedule GRHydroTransformADMToLocalBasis IN ADMBase_SetADMVars -{ - LANG: FORTRAN -} "Transform metric and shift to local tensor basis." + schedule GRHydroTransformADMToLocalBasis IN ADMBase_SetADMVars + { + LANG: C + } "Transform metric and shift to local tensor basis." -#schedule GRHydroTransformADMToLocalBasis IN CTG_Convert_to_ADM AFTER CTGBase_Convert_CTG_to_ADM -#{ -# LANG: FORTRAN -#} "Transform metric and shift to local tensor basis." + #schedule GRHydroTransformADMToLocalBasis IN CTG_Convert_to_ADM AFTER CTGBase_Convert_CTG_to_ADM + #{ + # LANG: C + #} "Transform metric and shift to local tensor basis." -#schedule GRHydroTransformADMToLocalBasis IN MoL_Step BEFORE MoL_CalcRHS -#{ -# LANG: FORTRAN -#} "Transform metric and shift to local tensor basis." + #schedule GRHydroTransformADMToLocalBasis IN MoL_Step BEFORE MoL_CalcRHS + #{ + # LANG: C + #} "Transform metric and shift to local tensor basis." -schedule GRHydroTransformPrimToGlobalBasis IN HydroBase_PostStep AFTER HydroBase_Con2Prim -{ - LANG: FORTRAN -} "Transform primitive vars to global tensor basis." + schedule GRHydroTransformPrimToGlobalBasis IN HydroBase_PostStep AFTER HydroBase_Con2Prim + { + LANG: C + } "Transform primitive vars to global tensor basis." +} schedule group GRHydroRHS IN HydroBase_RHS diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 index 0a1b3d8..a432d09 100644 --- a/src/GRHydro_Boundaries.F90 +++ b/src/GRHydro_Boundaries.F90 @@ -55,6 +55,7 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) integer, dimension(3) :: sym integer :: ierr = 0 integer :: itracer + CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates character(len=100) tracername character(len=100) tracerindex @@ -62,6 +63,8 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(1) = 1 sym(2) = 1 sym(3) = 1 + + general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::rho") call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::press") @@ -94,7 +97,11 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(2) = 1 sym(3) = 1 - call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[0]") + if (general_coordinates .ne. 0) then + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[0]") + else + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[0]") + endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[0]") if(evolve_mhd.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[0]") @@ -105,7 +112,11 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(2) = -1 sym(3) = 1 - call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[1]") + if (general_coordinates .ne. 0) then + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[1]") + else + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[1]") + endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[1]") if(evolve_mhd.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[1]") @@ -116,7 +127,11 @@ subroutine GRHydro_InitSymBound(CCTK_ARGUMENTS) sym(2) = 1 sym(3) = -1 - call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[2]") + if (general_coordinates .ne. 0) then + call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::lvel[2]") + else + call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::vel[2]") + endif call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::scon[2]") if(evolve_mhd.ne.0) then call SetCartSymVN(ierr, cctkGH, sym, "HydroBase::Bvec[2]") @@ -158,12 +173,15 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) integer :: ierr CCTK_REAL :: pi = 3.141569d0 integer :: i,j,k + CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates CCTK_INT, parameter :: faces=CCTK_ALL_FACES CCTK_INT, parameter :: ione=1 sw = GRHydro_stencil + general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) + !!$Flat boundaries if required if (CCTK_EQUALS(bound,"flat")) then @@ -181,11 +199,21 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) "HydroBase::press", "Flat") ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::eps", "Flat") - ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "Flat") - if(evolve_mhd.ne.0) then + if (general_coordinates .ne. 0) then ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") + "GRHydro::lvel", "Flat") + else + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", "Flat") + endif + if(evolve_mhd.ne.0) then + if (general_coordinates .ne. 0) then + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "Flat") + else + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "Flat") + endif ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "GRHydro::Bcons", "Flat") if(clean_divergence.ne.0) then @@ -236,8 +264,13 @@ subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & "HydroBase::vel", "None") if(evolve_mhd.ne.0) then - ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") + if (general_coordinates .ne. 0) then + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", "None") + else + ierr = Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", "None") + endif 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 3c51302..db03717 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -61,6 +61,14 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) CCTK_REAL :: local_min_tracer CCTK_REAL :: local_perc_ptol + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) @@ -68,6 +76,25 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) xpress = 0.0d0;xeps = 0.0d0;xtemp = 0.0d0;xye = 0.0d0 ! end EOS Omni vars + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") type2_bits = -1 @@ -101,11 +128,11 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) epsnegative = .false. - 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)) + det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ + g22(i,j,k),g23(i,j,k),g33(i,j,k)) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,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)) + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& + g23(i,j,k),g33(i,j,k)) !!$ if (det < 0.e0) then !!$ call CCTK_WARN(0, "nan produced (det negative)") @@ -135,7 +162,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 - lvel(i,j,k,:) = 0.d0 + vup(i,j,k,:) = 0.d0 w_lorentz(i,j,k) = 1.d0 if(evolve_temper.ne.0) then @@ -165,15 +192,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),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), & + scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2), & + vup(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),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),& + scon(i,j,k,2),scon(i,j,k,3),tau(i,j,k),Y_e_con(i,j,k),rho(i,j,k),vup(i,j,k,1),& + vup(i,j,k,2),vup(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 +213,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),& - lvel(i,j,k,1),lvel(i,j,k,2), lvel(i,j,k,3),eps(i,j,k),& + vup(i,j,k,1),vup(i,j,k,2), vup(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 +264,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 - lvel(i,j,k,:) = 0.d0 + vup(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 +286,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), 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), & + tau(i,j,k), rho(i,j,k), vup(i,j,k,1), vup(i,j,k,2), & + vup(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 @@ -1011,6 +1038,12 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) CCTK_REAL :: local_min_tracer + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) real*8 :: xpress,xtemp,xye,xeps @@ -1018,6 +1051,23 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0 ! end EOS Omni vars + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + end if + ! this is a poly call call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& GRHydro_rho_min,1.0d0,xtemp,xye,pmin,keyerr,anyerr) @@ -1048,18 +1098,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 * (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)) + gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) epsnegative = .false. @@ -1222,6 +1272,33 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) CCTK_REAL :: local_min_tracer ! character(len=400) :: warnline + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") @@ -1251,11 +1328,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(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 = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ + g22(i,j,k),g23(i,j,k),g33(i,j,k)) call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,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)) + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& + g23(i,j,k),g33(i,j,k)) if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -1277,8 +1354,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),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), & + scon(i,j,k,3),tau(i,j,k),rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2), & + vup(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)) @@ -1583,6 +1660,28 @@ subroutine Con2PrimBoundsPolytype(CCTK_ARGUMENTS) CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + end if call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere") @@ -1600,18 +1699,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 * (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)) + gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (g33(i,j,k) + g33(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) @@ -1677,6 +1776,28 @@ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS) uxxr, uxyr, uxzr, uyyr, uyzr, uzzr CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + end if nx = cctk_lsh(1) ny = cctk_lsh(2) @@ -1686,18 +1807,18 @@ subroutine Con2PrimBoundsTracer(CCTK_ARGUMENTS) do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - 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)) + gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (g33(i,j,k) + g33(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) @@ -1846,6 +1967,17 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS) integer :: i, j, k, nx, ny, nz character(len=100) warnline + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + vup => lvel + else + vup => vel + end if + ! call CCTK_INFO("Checking the C2P failure mask.") nx = cctk_lsh(1) @@ -1878,7 +2010,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: ',& - lvel(i,j,k,1),lvel(i,j,k,2),lvel(i,j,k,3) + vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3) call CCTK_WARN(0,warnline) !$OMP END CRITICAL diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90 index d1c022c..dda0070 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) 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(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(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) @@ -50,7 +50,15 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) CCTK_REAL :: det CCTK_REAL :: local_gamma - + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + !!$ Set up the fluid constants ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) @@ -64,6 +72,25 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr) + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => g33 + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + local_Gamma = 1.0d0 + xpress/xeps press = press_gf * poly_k_cgs * & (rho * inv_rho_gf)**local_Gamma @@ -78,10 +105,10 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ + g22(i,j,k),g23(i,j,k),g33(i,j,k)) + call prim2conpolytype(GRHydro_polytrope_handle,g11(i,j,k),g12(i,j,k),& + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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)) @@ -212,6 +239,15 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: Q character(len=100) infoline + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) real*8 :: xpress,xeps,xtemp,xye @@ -226,6 +262,25 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& 1.0d0,1.0d0,xtemp,xye,xpress,xeps,keyerr,anyerr) + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + local_Gamma = 1.0d0 + xpress/xeps local_K = xpress @@ -242,10 +297,10 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - 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 = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ + g22(i,j,k),g23(i,j,k),g33(i,j,k)) + call prim2con(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k),& + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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 ee6c0ad..474a2aa 100644 --- a/src/GRHydro_HLLE.F90 +++ b/src/GRHydro_HLLE.F90 @@ -50,6 +50,37 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) CCTK_INT :: type_bits, trivial, not_trivial + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + end if + if (flux_direction == 1) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & @@ -104,26 +135,26 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) !!$ left and right points. if (flux_direction == 1) then - avg_beta = 0.5d0 * (betaa(i+xoffset,j+yoffset,k+zoffset) + & - betaa(i,j,k)) + avg_beta = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + & + beta1(i,j,k)) else if (flux_direction == 2) then - avg_beta = 0.5d0 * (betab(i+xoffset,j+yoffset,k+zoffset) + & - betab(i,j,k)) + avg_beta = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + & + beta2(i,j,k)) else if (flux_direction == 3) then - avg_beta = 0.5d0 * (betac(i+xoffset,j+yoffset,k+zoffset) + & - betac(i,j,k)) + avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & + beta3(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 * (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)) + gxxh = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k)) + gxyh = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k)) + gxzh = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k)) + gyyh = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k)) + gyzh = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k)) + gzzh = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(i,j,k)) avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,\ gyyh,gyzh,gzzh) @@ -427,6 +458,37 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) CCTK_INT :: type_bits, trivial, not_trivial + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + end if + if (flux_direction == 1) then call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & @@ -467,26 +529,26 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) !!$ left and right points. if (flux_direction == 1) then - avg_beta = 0.5d0 * (betaa(i+xoffset,j+yoffset,k+zoffset) + & - betaa(i,j,k)) + avg_beta = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + & + beta1(i,j,k)) else if (flux_direction == 2) then - avg_beta = 0.5d0 * (betab(i+xoffset,j+yoffset,k+zoffset) + & - betab(i,j,k)) + avg_beta = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + & + beta2(i,j,k)) else if (flux_direction == 3) then - avg_beta = 0.5d0 * (betac(i+xoffset,j+yoffset,k+zoffset) + & - betac(i,j,k)) + avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & + beta3(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 * (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)) + gxxh = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k)) + gxyh = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k)) + gxzh = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k)) + gyyh = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k)) + gyzh = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k)) + gzzh = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(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 index 4dd6f47..0c26460 100644 --- a/src/GRHydro_Jacobian_state.c +++ b/src/GRHydro_Jacobian_state.c @@ -1,3 +1,5 @@ +#include <assert.h> + #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -15,13 +17,29 @@ GRHydro_check_Jacobian_state (CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS int ierr; + static CCTK_INT idxJacobian_state = -1, idxiJacobian_state = -1; + + if (!CCTK_IsImplementationActive("Coordinates")) + return; /* Multipatch infrastructure not active */ + + if (idxJacobian_state == -1) + { + idxJacobian_state = CCTK_VarIndex("Coordinates::jacobian_state"); + assert(idxJacobian_state >= 0); + } + + if (idxiJacobian_state == -1) + { + idxiJacobian_state = CCTK_VarIndex("Coordinates::inverse_jacobian_state"); + assert(idxiJacobian_state >= 0); + } - if (*jacobian_state == 0) + if (*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxJacobian_state) == 0) { CCTK_WARN(0, "GRHydro requires storage for Jacobians. Please tell your Coordinates implementation to store the Jacobians!"); } - if (*inverse_jacobian_state == 0) + if (*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian_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 a8058cc..2eecf6f 100644 --- a/src/GRHydro_PPMMReconstruct_drv.F90 +++ b/src/GRHydro_PPMMReconstruct_drv.F90 @@ -14,9 +14,9 @@ #include "SpaceMask.h" -#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(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(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) @@ -64,6 +64,41 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) CCTK_INT :: ierr + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + vup => vel + end if + allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with trivial_rp") @@ -152,8 +187,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& 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), betaa(:,j,k), alp(:,j,k),& + g11(:,j,k), g12(:,j,k), g13(:,j,k), g22(:,j,k), g23(:,j,k), & + g33(:,j,k), beta1(:,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, & @@ -167,8 +202,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) rhoplus(:,j,k),velxplus(:,j,k),velyplus(:,j,k),velzplus(:,j,k),& 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), betaa(:,j,k), alp(:,j,k),& + g11(:,j,k), g12(:,j,k), g13(:,j,k), g22(:,j,k), g23(:,j,k), & + g33(:,j,k), beta1(:,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, & @@ -223,8 +258,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& 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), betab(j,:,k), alp(j,:,k),& + g22(j,:,k), g23(j,:,k), g12(j,:,k), g33(j,:,k), g13(j,:,k), & + g11(j,:,k), beta2(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, & @@ -238,8 +273,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) rhoplus(j,:,k),velyplus(j,:,k),velzplus(j,:,k),velxplus(j,:,k),& 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), betab(j,:,k), alp(j,:,k),& + g22(j,:,k), g23(j,:,k), g12(j,:,k), g33(j,:,k), g13(j,:,k), & + g11(j,:,k), beta2(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, & @@ -293,8 +328,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& 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,:), betac(j,k,:), alp(j,k,:),& + g33(j,k,:), g13(j,k,:), g23(j,k,:), g11(j,k,:), g12(j,k,:), & + g22(j,k,:), beta3(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, & @@ -308,8 +343,8 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) rhoplus(j,k,:),velzplus(j,k,:),velxplus(j,k,:),velyplus(j,k,:),& 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,:), betac(j,k,:), alp(j,k,:),& + g33(j,k,:), g13(j,k,:), g23(j,k,:), g11(j,k,:), g12(j,k,:), & + g22(j,k,:), beta3(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 2ecbb71..aff5b7d 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) 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(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(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) @@ -55,6 +55,41 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) CCTK_INT :: ierr + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + vup => vel + end if + allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with trivial_rp") @@ -125,8 +160,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),& - gaa(:,j,k), gab(:,j,k), gac(:,j,k), gbb(:,j,k), gbc(:,j,k), & - gcc(:,j,k), betaa(:,j,k), alp(:,j,k),& + gxx(:,j,k), g12(:,j,k), g13(:,j,k), g22(:,j,k), g23(:,j,k), & + g33(:,j,k), beta1(:,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 +211,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),& - gbb(j,:,k), gbc(j,:,k), gab(j,:,k), gcc(j,:,k), gac(j,:,k), & - gaa(j,:,k), betab(j,:,k), alp(j,:,k),& + g22(j,:,k), g23(j,:,k), g12(j,:,k), g33(j,:,k), g13(j,:,k), & + gxx(j,:,k), beta2(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 +262,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,:),& - gcc(j,k,:), gac(j,k,:), gbc(j,k,:), gaa(j,k,:), gab(j,k,:), & - gbb(j,k,:), betac(j,k,:), alp(j,k,:),& + g33(j,k,:), g13(j,k,:), g23(j,k,:), gxx(j,k,:), g12(j,k,:), & + g22(j,k,:), beta3(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_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90 index 181a1cb..9310448 100644 --- a/src/GRHydro_ParamCheck.F90 +++ b/src/GRHydro_ParamCheck.F90 @@ -35,6 +35,8 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + integer :: coordinates_is_active + if (GRHydro_stencil > minval(cctk_nghostzones)) then call CCTK_PARAMWARN("The stencil is larger than the number of ghost zones. Answer will be dependent on processor number...") end if @@ -104,6 +106,13 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) evolve_temper = 0 endif + call CCTK_IsImplementationActive(coordinates_is_active, "Coordinates") + ! this test is somewhat overzealous but we cannot access + ! coordinates::general_coordinates yet since it is only set in BaseGrid + if (CCTK_Equals(riemann_solver,"HLLE").eq.0.and.coordinates_is_active.ne.0) then + call CCTK_WARN(0, "There is currently no Riemann solver other than HLLE compatible with multipatch!") + end if + if (CCTK_EQUALS(riemann_solver,"Roe").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then call CCTK_PARAMWARN("Roe solver is not implemented yet for MHD") end if diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 2bd2d29..ee82403 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) 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(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(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,70 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& - gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_detr CCTK_REAL :: xtemp(1) + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr,& - !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, & - !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr) + !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & + !$OMP g11r,g12r,g13r,g22r,g23r,g33r) 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 - 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)) + g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + g33r = 0.5d0 * (g33(i,j,k) + g33(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) + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) - call prim2con(GRHydro_eos_handle, gaal,gabl,gacl,gbbl,& - gbcl,gccl, & + call prim2con(GRHydro_eos_handle, g11l,g12l,g13l,g22l,& + g23l,g33l, & 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, gaar,gabr,gacr,& - gbbr,gbcr,gccr, & + call prim2con(GRHydro_eos_handle, g11r,g12r,g13r,& + g22r,g23r,g33r, & 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 +120,27 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) !$OMP END PARALLEL DO else !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& - !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, & - !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr) + !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & + !$OMP g11r,g12r,g13r,g22r,g23r,g33r) 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 - 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)) + g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + g33r = 0.5d0 * (g33(i,j,k) + g33(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) + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -124,8 +150,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),gaal,gabl,gacl,gbbl,& - gbcl,gccl, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11l,g12l,g13l,g22l,& + g23l,g33l, & 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 +167,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),gaar,gabr,gacr,& - gbbr,gbcr,gccr, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11r,g12r,g13r,& + g22r,g23r,g33r, & 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),& @@ -347,18 +373,44 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) CCTK_REAL :: det character(len=512) :: warnline + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, k, det) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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)) + det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ + g22(i,j,k),g23(i,j,k),g33(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),& + call prim2con(GRHydro_eos_handle,g11(i,j,k),& + g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(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 +424,13 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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)) + det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ + g22(i,j,k),g23(i,j,k),g33(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),gaa(i,j,k),& - gab(i,j,k),gac(i,j,k),& - gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),& + x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),& + g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(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 +481,66 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& - gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if - !$OMP PARALLEL DO PRIVATE(i, j, k, gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& - !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r,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 - 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, & + g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + + call prim2conpolytype(GRHydro_eos_handle, g11l,g12l,g13l,& + g22l,g23l,g33l, & 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, gaar,gabr,gacr,& - gbbr,gbcr,gccr, & + call prim2conpolytype(GRHydro_eos_handle, g11r,g12r,g13r,& + g22r,g23r,g33r, & 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),& @@ -580,17 +658,43 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_REAL :: det + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if !$OMP PARALLEL DO PRIVATE(i, j, k, det) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - 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)) + det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ + g22(i,j,k),g23(i,j,k),g33(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),& + call prim2conpolytype(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k),& + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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 +741,68 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& - gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + end if 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 - 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) + g11l = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + g12l = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + g13l = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + g22l = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + g23l = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + g33l = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + g11r = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + g12r = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + g13r = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + g22r = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) + avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * & sqrt(avg_detr) * rhoplus(i,j,k) / & sqrt(1.d0 - & - (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) ) ) ) + (g11r * velxplus(i,j,k)**2 + & + g22r * velyplus(i,j,k)**2 + & + g33r * velzplus(i,j,k)**2 + & + 2.d0 * (g12r * velxplus(i,j,k) * velyplus(i,j,k) + & + g13r * velxplus(i,j,k) * velzplus(i,j,k) + & + g23r * 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 - & - (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) ) ) ) + (g11l * velxminus(i,j,k)**2 + & + g22l * velyminus(i,j,k)**2 + & + g33l * velzminus(i,j,k)**2 + & + 2.d0 * (g12l * velxminus(i,j,k) * velyminus(i,j,k) + & + g13l * velxminus(i,j,k) * velzminus(i,j,k) + & + g23l * 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 5e91721..43346ed 100644 --- a/src/GRHydro_RegisterVars.cc +++ b/src/GRHydro_RegisterVars.cc @@ -13,6 +13,8 @@ using namespace std; +extern "C" CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); + // Utility functions to register variables with MoL // Note: We could check for the return value here, but MoL issues a // level 0 warning in that case anyway. If that changes in the @@ -37,6 +39,8 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; + const int general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH); + // We need some aliased functions, so we first check if they are available string needed_funs[5] = {"MoLRegisterEvolvedGroup", "MoLRegisterConstrainedGroup", @@ -55,11 +59,13 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) register_constrained("HydroBase::eps"); register_constrained("HydroBase::w_lorentz"); register_constrained("HydroBase::vel"); - register_constrained("GRHydro::lvel"); + if (general_coordinates) { + register_constrained("GRHydro::lvel"); - register_constrained("grhydro::local_shift"); - register_constrained("grhydro::local_metric"); - register_constrained("grhydro::local_extrinsic_curvature"); + register_constrained("grhydro::local_shift"); + register_constrained("grhydro::local_metric"); + register_constrained("grhydro::local_extrinsic_curvature"); + } if (CCTK_EQUALS(evolution_method, "GRHydro")) { @@ -69,7 +75,9 @@ extern "C"void GRHydro_Register(CCTK_ARGUMENTS) if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro")) { register_constrained("HydroBase::Bvec"); - register_constrained("GRHydro::lBvec"); + if (general_coordinates) { + 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 969e88f..71dd75a 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -27,9 +27,9 @@ #include "cctk_Arguments.h" #include "GRHydro_Macros.h" -#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) +#define vela(i,j,k) vup(i,j,k,1) +#define velb(i,j,k) vup(i,j,k,2) +#define velc(i,j,k) vup(i,j,k,3) /*@@ @routine SourceTerms @date Sat Jan 26 02:04:21 2002 @@ -76,6 +76,41 @@ subroutine SourceTerms(CCTK_ARGUMENTS) logical, allocatable, dimension (:,:,:) :: force_spatial_second_order + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + vup => vel + end if + one = 1.0d0 two = 2.0d0 half = 0.5d0 @@ -156,12 +191,12 @@ subroutine SourceTerms(CCTK_ARGUMENTS) !!$ Set the metric terms. - 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) + localgaa = g11(i,j,k) + localgab = g12(i,j,k) + localgac = g13(i,j,k) + localgbb = g22(i,j,k) + localgbc = g23(i,j,k) + localgcc = g33(i,j,k) det = SPATIAL_DETERMINANT(localgaa, localgab, localgac,\ localgbb, localgbc, localgcc) @@ -175,37 +210,37 @@ subroutine SourceTerms(CCTK_ARGUMENTS) rhoenthalpyW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k))*& w_lorentz(i,j,k)**2 - shifta = betaa(i,j,k) - shiftb = betab(i,j,k) - shiftc = betac(i,j,k) + shifta = beta1(i,j,k) + shiftb = beta2(i,j,k) + shiftc = beta3(i,j,k) if (local_spatial_order .eq. 2) then - da_betaa = DIFF_X_2(betaa) - da_betab = DIFF_X_2(betab) - da_betac = DIFF_X_2(betac) + da_betaa = DIFF_X_2(beta1) + da_betab = DIFF_X_2(beta2) + da_betac = DIFF_X_2(beta3) - db_betaa = DIFF_Y_2(betaa) - db_betab = DIFF_Y_2(betab) - db_betac = DIFF_Y_2(betac) + db_betaa = DIFF_Y_2(beta1) + db_betab = DIFF_Y_2(beta2) + db_betac = DIFF_Y_2(beta3) - dc_betaa = DIFF_Z_2(betaa) - dc_betab = DIFF_Z_2(betab) - dc_betac = DIFF_Z_2(betac) + dc_betaa = DIFF_Z_2(beta1) + dc_betab = DIFF_Z_2(beta2) + dc_betac = DIFF_Z_2(beta3) else - da_betaa = DIFF_X_4(betaa) - da_betab = DIFF_X_4(betab) - da_betac = DIFF_X_4(betac) + da_betaa = DIFF_X_4(beta1) + da_betab = DIFF_X_4(beta2) + da_betac = DIFF_X_4(beta3) - db_betaa = DIFF_Y_4(betaa) - db_betab = DIFF_Y_4(betab) - db_betac = DIFF_Y_4(betac) + db_betaa = DIFF_Y_4(beta1) + db_betab = DIFF_Y_4(beta2) + db_betac = DIFF_Y_4(beta3) - dc_betaa = DIFF_Z_4(betaa) - dc_betab = DIFF_Z_4(betab) - dc_betac = DIFF_Z_4(betac) + dc_betaa = DIFF_Z_4(beta1) + dc_betab = DIFF_Z_4(beta2) + dc_betac = DIFF_Z_4(beta3) end if @@ -261,45 +296,45 @@ subroutine SourceTerms(CCTK_ARGUMENTS) if (local_spatial_order .eq. 2) then - 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) + da_gaa = DIFF_X_2(g11) + da_gab = DIFF_X_2(g12) + da_gac = DIFF_X_2(g13) + da_gbb = DIFF_X_2(g22) + da_gbc = DIFF_X_2(g23) + da_gcc = DIFF_X_2(g33) + db_gaa = DIFF_Y_2(g11) + db_gab = DIFF_Y_2(g12) + db_gac = DIFF_Y_2(g13) + db_gbb = DIFF_Y_2(g22) + db_gbc = DIFF_Y_2(g23) + db_gcc = DIFF_Y_2(g33) + dc_gaa = DIFF_Z_2(g11) + dc_gab = DIFF_Z_2(g12) + dc_gac = DIFF_Z_2(g13) + dc_gbb = DIFF_Z_2(g22) + dc_gbc = DIFF_Z_2(g23) + dc_gcc = DIFF_Z_2(g33) else - 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) + da_gaa = DIFF_X_4(g11) + da_gab = DIFF_X_4(g12) + da_gac = DIFF_X_4(g13) + da_gbb = DIFF_X_4(g22) + da_gbc = DIFF_X_4(g23) + da_gcc = DIFF_X_4(g33) + db_gaa = DIFF_Y_4(g11) + db_gab = DIFF_Y_4(g12) + db_gac = DIFF_Y_4(g13) + db_gbb = DIFF_Y_4(g22) + db_gbc = DIFF_Y_4(g23) + db_gcc = DIFF_Y_4(g33) + dc_gaa = DIFF_Z_4(g11) + dc_gab = DIFF_Z_4(g12) + dc_gac = DIFF_Z_4(g13) + dc_gbb = DIFF_Z_4(g22) + dc_gbc = DIFF_Z_4(g23) + dc_gcc = DIFF_Z_4(g33) end if diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90 index 96a71db..97db410 100644 --- a/src/GRHydro_TVDReconstruct_drv.F90 +++ b/src/GRHydro_TVDReconstruct_drv.F90 @@ -14,9 +14,9 @@ #include "SpaceMask.h" -#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(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(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) @@ -62,6 +62,17 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) CCTK_INT :: ierr + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET vel, lvel + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + vup => lvel + else + vup => vel + end if + allocate(trivial_rp(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),STAT=ierr) if (ierr .ne. 0) then call CCTK_WARN(0, "Allocation problems with trivial_rp") @@ -155,11 +166,11 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) rho, rhoplus, rhominus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - lvel(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) + vup(:,:,:,1), velxplus, velxminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - lvel(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) + vup(:,:,:,2), velyplus, velyminus, trivial_rp, hydro_excision_mask) call tvdreconstruct(nx, ny, nz, xoffset, yoffset, zoffset, & - lvel(:,:,:,3), velzplus, velzminus, trivial_rp, hydro_excision_mask) + vup(:,:,:,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 diff --git a/src/GRHydro_TransformTensorBasis.F90 b/src/GRHydro_TransformTensorBasis.F90 deleted file mode 100644 index 1aeda46..0000000 --- a/src/GRHydro_TransformTensorBasis.F90 +++ /dev/null @@ -1,211 +0,0 @@ -#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_TransformTensorBasis.c b/src/GRHydro_TransformTensorBasis.c new file mode 100644 index 0000000..f36cb36 --- /dev/null +++ b/src/GRHydro_TransformTensorBasis.c @@ -0,0 +1,438 @@ +#include <assert.h> + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#define SQR(x) ((x)*(x)) + +/* helper routine to let Fortran code test for the presence of Multipatch. + * Required since Fortran cannot access grid scalars via pointers returned by CCTK_VarDataPtrI + * Returns 0 when no general coordinates are used + * Returns 1 when general coordinates are used + */ +CCTK_FCALL CCTK_INT CCTK_FNAME(GRHydro_UseGeneralCoordinates)(CCTK_POINTER ptr_to_cctkGH); +CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); + +CCTK_FCALL CCTK_INT CCTK_FNAME(GRHydro_UseGeneralCoordinates)(CCTK_POINTER ptr_to_cctkGH) +{ + return GRHydro_UseGeneralCoordinates(*(cGH **)ptr_to_cctkGH); +} + +CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH) +{ + static CCTK_INT idxGeneral_coordinates = -1; + static CCTK_INT coordinates_active_set = 0; + static CCTK_INT coordinates_active; + + if (!coordinates_active_set) { + coordinates_active = CCTK_IsImplementationActive("Coordinates"); + coordinates_active_set = 1; + } + + if (idxGeneral_coordinates == -1) + { + idxGeneral_coordinates = CCTK_VarIndex("Coordinates::general_coordinates"); + assert(!coordinates_active || idxGeneral_coordinates >= 0); + } + + /* are coordinates relly not Cartesian? */ + return coordinates_active && *(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxGeneral_coordinates); +} + +void GRHydroTransformADMToLocalBasis(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int i,j,k; + static CCTK_INT idxJacobian = -1, idxiJacobian = -1; + static CCTK_INT idxJacobian_state = -1, idxiJacobian_state = -1; + static CCTK_INT idxGeneral_coordinates = -1; + CCTK_REAL *J11, *J12, *J13, *J21, *J22, *J23, *J31, *J32, *J33; + CCTK_REAL *iJ11, *iJ12, *iJ13, *iJ21, *iJ22, *iJ23, *iJ31, *iJ32, *iJ33; + + if (idxGeneral_coordinates == -1) + { + idxGeneral_coordinates = CCTK_VarIndex("Coordinates::general_coordinates"); + assert(idxGeneral_coordinates >= 0); + } + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxGeneral_coordinates)) + return; /* corodinates are Cartesian */ + + if (idxJacobian == -1) + { + idxJacobian = CCTK_FirstVarIndex("Coordinates::jacobian"); + assert(idxJacobian >= 0); + } + + if (idxiJacobian == -1) + { + idxiJacobian = CCTK_FirstVarIndex("Coordinates::inverse_jacobian"); + assert(idxiJacobian >= 0); + } + + if (idxJacobian_state == -1) + { + idxJacobian_state = CCTK_VarIndex("Coordinates::jacobian_state"); + assert(idxJacobian_state >= 0); + } + + if (idxiJacobian_state == -1) + { + idxiJacobian_state = CCTK_VarIndex("Coordinates::inverse_jacobian_state"); + assert(idxiJacobian_state >= 0); + } + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxJacobian_state)) + CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!"); + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian_state)) + CCTK_WARN(0,"No storage for inverse Jacobians allocated! Tell your coordinates implemetation to store inverse Jacobians!"); + + /* Cactus guarantess that variables in a group are consecutive, and + * ReflectionSymmetry etc. require the order we use */ + J11 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+0); + J12 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+1); + J13 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+2); + J21 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+3); + J22 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+4); + J23 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+5); + J31 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+6); + J32 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+7); + J33 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+8); + assert(J11 && J12 && J13 && J21 && J22 && J23 && J31 && J32 && J33); + iJ11 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+0); + iJ12 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+1); + iJ13 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+2); + iJ21 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+3); + iJ22 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+4); + iJ23 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+5); + iJ31 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+6); + iJ32 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+7); + iJ33 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+8); + assert(iJ11 && iJ12 && iJ13 && iJ21 && iJ22 && iJ23 && iJ31 && iJ32 && iJ33); + + #pragma omp parallel for private (i,j,k) + for (k=0 ; k<cctk_lsh[2] ; k++) + { + for (j=0 ; j<cctk_lsh[1] ; j++) + { + for (i=0 ; i<cctk_lsh[0] ; i++) + { + CCTK_INT idx = CCTK_GFINDEX3D(cctkGH, i,j,k); + + gaa[idx] = SQR(iJ11[idx]) * gxx[idx] + + 2.0 * iJ11[idx] * iJ21[idx] * gxy[idx] + + 2.0 * iJ11[idx] * iJ31[idx] * gxz[idx] + SQR(iJ21[idx]) * gyy[idx] + + 2.0 * iJ21[idx] * iJ31[idx] * gyz[idx] + SQR(iJ31[idx]) * gzz[idx]; + + gab[idx] = iJ11[idx] * iJ12[idx] * gxx[idx] + + (iJ11[idx] * iJ22[idx] + iJ21[idx] * iJ12[idx]) * gxy[idx] + + (iJ11[idx] * iJ32[idx] + iJ31[idx] * iJ12[idx]) * gxz[idx] + + iJ21[idx] * iJ22[idx] * gyy[idx] + (iJ21[idx] * iJ32[idx] + + iJ31[idx] * iJ22[idx]) * gyz[idx] + iJ31[idx] * iJ32[idx] * gzz[idx]; + + gac[idx] = iJ11[idx] * iJ13[idx] * gxx[idx] + (iJ11[idx] * iJ23[idx] + + iJ21[idx] * iJ13[idx]) * gxy[idx] + + (iJ11[idx] * iJ33[idx] + iJ31[idx] * iJ13[idx]) * gxz[idx] + + iJ21[idx] * iJ23[idx] * gyy[idx] + + (iJ21[idx] * iJ33[idx] + iJ31[idx] * iJ23[idx]) * gyz[idx] + + iJ31[idx] * iJ33[idx] * gzz[idx]; + + gbb[idx] = SQR(iJ12[idx]) * gxx[idx] + + 2.0 * iJ12[idx] * iJ22[idx] * gxy[idx] + + 2.0 * iJ12[idx] * iJ32[idx] * gxz[idx] + SQR(iJ22[idx]) * gyy[idx] + + 2.0 * iJ22[idx] * iJ32[idx] * gyz[idx] + SQR(iJ32[idx]) * gzz[idx]; + + gbc[idx] = iJ12[idx] * iJ13[idx] * gxx[idx] + + (iJ12[idx] * iJ23[idx] + iJ22[idx] * iJ13[idx]) * gxy[idx] + + (iJ12[idx] * iJ33[idx] + iJ32[idx] * iJ13[idx]) * gxz[idx] + + iJ22[idx] * iJ23[idx] * gyy[idx] + + (iJ22[idx]* iJ33[idx] + iJ32[idx] * iJ23[idx]) * gyz[idx] + + iJ32[idx] * iJ33[idx] * gzz[idx]; + + gcc[idx] = SQR(iJ13[idx]) * gxx[idx] + + 2.0 * iJ13[idx] * iJ23[idx] * gxy[idx] + + 2.0 * iJ13[idx] * iJ33[idx] * gxz[idx] + + SQR(iJ23[idx]) * gyy[idx] + + 2.0 * iJ23[idx] * iJ33[idx] * gyz[idx] + + SQR(iJ33[idx]) * gzz[idx]; + + /* Transform extrinsic curvature from global to local basis. + * Since extrinsic has covariant indices, use inverse Jacobian */ + kaa[idx] = SQR(iJ11[idx]) * kxx[idx] + + 2.0 * iJ11[idx] * iJ21[idx] * kxy[idx] + + 2.0 * iJ11[idx] * iJ31[idx] * kxz[idx] + + SQR(iJ21[idx]) * kyy[idx] + + 2.0 * iJ21[idx] * iJ31[idx] * kyz[idx] + + SQR(iJ31[idx]) * kzz[idx]; + + kab[idx] = iJ11[idx] * iJ12[idx] * kxx[idx] + + (iJ11[idx] * iJ22[idx] + iJ21[idx] * iJ12[idx]) * kxy[idx] + + (iJ11[idx] * iJ32[idx] + iJ31[idx] * iJ12[idx]) * kxz[idx] + + iJ21[idx] * iJ22[idx] * kyy[idx] + + (iJ21[idx] * iJ32[idx] + iJ31[idx] * iJ22[idx]) * kyz[idx] + + iJ31[idx] * iJ32[idx] * kzz[idx]; + + + kac[idx] = iJ11[idx] * iJ13[idx] * kxx[idx] + + (iJ11[idx] * iJ23[idx] + iJ21[idx] * iJ13[idx]) * kxy[idx] + + (iJ11[idx] * iJ33[idx] + iJ31[idx] * iJ13[idx]) * kxz[idx] + + iJ21[idx] * iJ23[idx] * kyy[idx] + + (iJ21[idx] * iJ33[idx] + iJ31[idx] * iJ23[idx]) * kyz[idx] + + iJ31[idx] * iJ33[idx] * kzz[idx]; + + kbb[idx] = SQR(iJ12[idx]) * kxx[idx] + + 2.0 * iJ12[idx] * iJ22[idx] * kxy[idx] + + 2.0 * iJ12[idx] * iJ32[idx] * kxz[idx] + + SQR(iJ22[idx]) * kyy[idx] + + 2.0 * iJ22[idx] * iJ32[idx] * kyz[idx] + + SQR(iJ32[idx]) * kzz[idx]; + + kbc[idx] = iJ12[idx] * iJ13[idx] * kxx[idx] + + (iJ12[idx] * iJ23[idx] + iJ22[idx] * iJ13[idx]) * kxy[idx] + + (iJ12[idx] * iJ33[idx] + iJ32[idx] * iJ13[idx]) * kxz[idx] + + iJ22[idx] * iJ23[idx] * kyy[idx] + + (iJ22[idx] * iJ33[idx] + iJ32[idx] * iJ23[idx]) * kyz[idx] + + iJ32[idx] * iJ33[idx] * kzz[idx]; + + kcc[idx] = SQR(iJ13[idx]) * kxx[idx] + + 2.0 * iJ13[idx] * iJ23[idx] * kxy[idx] + + 2.0 * iJ13[idx] * iJ33[idx] * kxz[idx] + + SQR(iJ23[idx]) * kyy[idx] + + 2.0 * iJ23[idx] * iJ33[idx] * kyz[idx] + + SQR(iJ33[idx]) * kzz[idx]; + + /* Transform shift from global to local basis. + * Since shift has contravariant index, use Jacobian */ + betaa[idx] = betax[idx]*J11[idx] + betay[idx]*J12[idx] + betaz[idx]*J13[idx]; + betab[idx] = betax[idx]*J21[idx] + betay[idx]*J22[idx] + betaz[idx]*J23[idx]; + betac[idx] = betax[idx]*J31[idx] + betay[idx]*J32[idx] + betaz[idx]*J33[idx]; + } + } + } + +} + + + +void GRHydroTransformPrimToLocalBasis(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT i,j,k; + static CCTK_INT idxJacobian = -1, idxiJacobian = -1; + static CCTK_INT idxGeneral_coordinates = -1; + static CCTK_INT idxJacobian_state = -1, idxiJacobian_state = -1; + CCTK_REAL *J11, *J12, *J13, *J21, *J22, *J23, *J31, *J32, *J33; + CCTK_REAL *iJ11, *iJ12, *iJ13, *iJ21, *iJ22, *iJ23, *iJ31, *iJ32, *iJ33; + + if (!CCTK_IsImplementationActive("Coordinates")) + return; /* Multipatch infrastructure not active */ + + if (idxGeneral_coordinates == -1) + { + idxGeneral_coordinates = CCTK_VarIndex("Coordinates::general_coordinates"); + assert(idxGeneral_coordinates >= 0); + } + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxGeneral_coordinates)) + return; /* corodinates are Cartesian */ + + if (idxJacobian == -1) + { + idxJacobian = CCTK_FirstVarIndex("Coordinates::jacobian"); + assert(idxJacobian >= 0); + } + + if (idxiJacobian == -1) + { + idxiJacobian = CCTK_FirstVarIndex("Coordinates::inverse_jacobian"); + assert(idxiJacobian >= 0); + } + + if (idxJacobian_state == -1) + { + idxJacobian_state = CCTK_VarIndex("Coordinates::jacobian_state"); + assert(idxJacobian_state >= 0); + } + + if (idxiJacobian_state == -1) + { + idxiJacobian_state = CCTK_VarIndex("Coordinates::inverse_jacobian_state"); + assert(idxiJacobian_state >= 0); + } + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxJacobian_state)) + CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!"); + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian_state)) + CCTK_WARN(0,"No storage for inverse Jacobians allocated! Tell your coordinates implemetation to store inverse Jacobians!"); + + /* CCTK_INFO("Transforming primitives to local basis."); */ + + /* Cactus guarantess that variables in a group are consecutive, and + * ReflectionSymmetry etc. require the order we use */ + J11 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+0); + J12 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+1); + J13 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+2); + J21 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+3); + J22 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+4); + J23 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+5); + J31 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+6); + J32 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+7); + J33 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+8); + assert(J11 && J12 && J13 && J21 && J22 && J23 && J31 && J32 && J33); + iJ11 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+0); + iJ12 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+1); + iJ13 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+2); + iJ21 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+3); + iJ22 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+4); + iJ23 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+5); + iJ31 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+6); + iJ32 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+7); + iJ33 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+8); + assert(iJ11 && iJ12 && iJ13 && iJ21 && iJ22 && iJ23 && iJ31 && iJ32 && iJ33); + + #pragma omp parallel for private (i,j,k) + for (k=0 ; k<cctk_lsh[2] ; k++) + { + for (j=0 ; j<cctk_lsh[1] ; j++) + { + for (i=0 ; i<cctk_lsh[0] ; i++) + { + CCTK_INT idx = CCTK_GFINDEX3D(cctkGH, i,j,k); + CCTK_INT idx1 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,0); + CCTK_INT idx2 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,1); + CCTK_INT idx3 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,2); + + /* Transform primitive velocity from global to local basis. + * Since velocity has contravariant index, use Jacobian */ + lvel[idx1] = vel[idx1]*J11[idx] + vel[idx2]*J12[idx] + vel[idx3]*J13[idx]; + lvel[idx2] = vel[idx1]*J21[idx] + vel[idx2]*J22[idx] + vel[idx3]*J23[idx]; + lvel[idx3] = vel[idx1]*J31[idx] + vel[idx2]*J32[idx] + vel[idx3]*J33[idx]; + + /* Transform primitive B-field from global to local basis. + * Since B-field has contravariant index, use Jacobian */ +/* lBvec[idx1] = Bvec[idx1]*J11[idx] + Bvec[idx2]*J12[idx] + Bvec[idx3]*J13[idx]; */ +/* lBvec[idx2] = Bvec[idx1]*J21[idx] + Bvec[idx2]*J22[idx] + Bvec[idx3]*J23[idx]; */ +/* lBvec[idx3] = Bvec[idx1]*J31[idx] + Bvec[idx2]*J32[idx] + Bvec[idx3]*J33[idx]; */ + } + } + } + +} + + + +void GRHydroTransformPrimToGlobalBasis(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT i,j,k; + + static CCTK_INT idxJacobian = -1, idxiJacobian = -1; + static CCTK_INT idxJacobian_state = -1, idxiJacobian_state = -1; + static CCTK_INT idxGeneral_coordinates = -1; + CCTK_REAL *J11, *J12, *J13, *J21, *J22, *J23, *J31, *J32, *J33; + CCTK_REAL *iJ11, *iJ12, *iJ13, *iJ21, *iJ22, *iJ23, *iJ31, *iJ32, *iJ33; + + if (idxGeneral_coordinates == -1) + { + idxGeneral_coordinates = CCTK_VarIndex("Coordinates::general_coordinates"); + assert(idxGeneral_coordinates >= 0); + } + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxGeneral_coordinates)) + return; /* corodinates are Cartesian */ + + if (idxJacobian == -1) + { + idxJacobian = CCTK_FirstVarIndex("Coordinates::jacobian"); + assert(idxJacobian >= 0); + } + + if (idxiJacobian == -1) + { + idxiJacobian = CCTK_FirstVarIndex("Coordinates::inverse_jacobian"); + assert(idxiJacobian >= 0); + } + + if (idxJacobian_state == -1) + { + idxJacobian_state = CCTK_VarIndex("Coordinates::jacobian_state"); + assert(idxJacobian_state >= 0); + } + + if (idxiJacobian_state == -1) + { + idxiJacobian_state = CCTK_VarIndex("Coordinates::inverse_jacobian_state"); + assert(idxiJacobian_state >= 0); + } + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxJacobian_state)) + CCTK_WARN(0,"No storage for Jacobians allocated! Tell your coordinates implemetation to store Jacobians!"); + + if (!*(CCTK_INT *)CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian_state)) + CCTK_WARN(0,"No storage for inverse Jacobians allocated! Tell your coordinates implemetation to store inverse Jacobians!"); + + /* CCTK_INFO("Transforming primitives to global basis."); */ + + /* Cactus guarantess that variables in a group are consecutive, and + * ReflectionSymmetry etc. require the order we use */ + J11 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+0); + J12 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+1); + J13 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+2); + J21 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+3); + J22 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+4); + J23 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+5); + J31 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+6); + J32 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+7); + J33 = CCTK_VarDataPtrI(cctkGH, 0, idxJacobian+8); + assert(J11 && J12 && J13 && J21 && J22 && J23 && J31 && J32 && J33); + iJ11 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+0); + iJ12 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+1); + iJ13 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+2); + iJ21 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+3); + iJ22 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+4); + iJ23 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+5); + iJ31 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+6); + iJ32 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+7); + iJ33 = CCTK_VarDataPtrI(cctkGH, 0, idxiJacobian+8); + assert(iJ11 && iJ12 && iJ13 && iJ21 && iJ22 && iJ23 && iJ31 && iJ32 && iJ33); + + #pragma omp parallel for private (i,j,k) + for (k=0 ; k<cctk_lsh[2] ; k++) + { + for (j=0 ; j<cctk_lsh[1] ; j++) + { + for (i=0 ; i<cctk_lsh[0] ; i++) + { + CCTK_INT idx = CCTK_GFINDEX3D(cctkGH, i,j,k); + CCTK_INT idx1 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,0); + CCTK_INT idx2 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,1); + CCTK_INT idx3 = CCTK_VECTGFINDEX3D(cctkGH, i,j,k,2); + + /* Transform primitive velocity from local to global basis. + * Since velocity has contravariant index, use inverse Jacobian */ + vel[idx1] = lvel[idx1]*iJ11[idx] + lvel[idx2]*iJ12[idx] + lvel[idx3]*iJ13[idx]; + vel[idx2] = lvel[idx1]*iJ21[idx] + lvel[idx2]*iJ22[idx] + lvel[idx3]*iJ23[idx]; + vel[idx3] = lvel[idx1]*iJ31[idx] + lvel[idx2]*iJ32[idx] + lvel[idx3]*iJ33[idx]; + + } + } + } + +} + + + + + diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90 index 0bb092c..c8a829e 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) 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) +#define velx(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(i,j,k,3) +#define velx_p(i,j,k) vup_p(i,j,k,1) +#define vely_p(i,j,k) vup_p(i,j,k,2) +#define velz_p(i,j,k) vup_p(i,j,k,3) +#define velx_p_p(i,j,k) vup_p_p(i,j,k,1) +#define vely_p_p(i,j,k) vup_p_p(i,j,k,2) +#define velz_p_p(i,j,k) vup_p_p(i,j,k,3) subroutine GRHydroUpdateAtmosphereMask(CCTK_ARGUMENTS) @@ -220,6 +220,15 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) CCTK_REAL :: det, psi4pt CCTK_INT :: type_bits, atmosphere, not_atmosphere + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) real*8 :: xpress,xeps,xtemp,xye @@ -227,6 +236,25 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0 ! end EOS Omni vars + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",& "in_atmosphere") @@ -246,8 +274,8 @@ 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(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 = SPATIAL_DETERMINANT(g11(i,j,k), g12(i,j,k), g13(i,j,k), \ + g22(i,j,k), g23(i,j,k), g33(i,j,k)) if(evolve_temper.ne.0) then ! ! set the temperature to be relatively low @@ -259,8 +287,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),gaa(i,j,k),gab(i,j,k),& - gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),g12(i,j,k),& + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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,8 +296,8 @@ 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, & - 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, & + g11(i,j,k), g12(i,j,k), g13(i,j,k), & + g22(i,j,k), g23(i,j,k), g33(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)) @@ -304,6 +332,29 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) CCTK_INT :: type_bits, atmosphere, not_atmosphere CCTK_INT :: eos_handle + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + TARGET gaa_p, gab_p, gac_p, gbb_p, gbc_p, gcc_p + TARGET gxx_p, gxy_p, gxz_p, gyy_p, gyz_p, gzz_p + TARGET lvel_p, vel_p + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11_p => Null(), g12_p => Null(), & + g13_p => Null(), g22_p => Null(), & + g23_p => Null(), g33_p => Null() + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup_p => Null() + TARGET gaa_p_p, gab_p_p, gac_p_p, gbb_p_p, gbc_p_p, gcc_p_p + TARGET gxx_p_p, gxy_p_p, gxz_p_p, gyy_p_p, gyz_p_p, gzz_p_p + TARGET lvel_p_p, vel_p_p + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11_p_p => Null(), g12_p_p => Null(), & + g13_p_p => Null(), g22_p_p => Null(), & + g23_p_p => Null(), g33_p_p => Null() + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup_p_p => Null() + ! begin EOS Omni vars integer :: n,keytemp,anyerr,keyerr(1) real*8 :: xpress,xeps,xtemp,xye @@ -311,6 +362,64 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0 ! end EOS Omni vars + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + end if + if (timelevels .gt. 1) then + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11_p => gaa_p + g12_p => gab_p + g13_p => gac_p + g22_p => gbb_p + g23_p => gbc_p + g33_p => gcc_p + vup_p => lvel_p + else + g11_p => gxx_p + g12_p => gxy_p + g13_p => gxz_p + g22_p => gyy_p + g23_p => gyz_p + g33_p => gzz_p + vup_p => vel_p + end if + end if + if (timelevels .gt. 2) then + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11_p_p => gaa_p_p + g12_p_p => gab_p_p + g13_p_p => gac_p_p + g22_p_p => gbb_p_p + g23_p_p => gbc_p_p + g33_p_p => gcc_p_p + vup_p_p => lvel_p_p + else + g11_p_p => gxx_p_p + g12_p_p => gxy_p_p + g13_p_p => gxz_p_p + g22_p_p => gyy_p_p + g23_p_p => gyz_p_p + g33_p_p => gzz_p_p + vup_p_p => vel_p_p + end if + end if + + eos_handle = GRHydro_polytrope_handle call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") @@ -330,8 +439,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) if (rho(i,j,k) .le. rho_min) then rho(i,j,k) = rho_min - 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)) + det = SPATIAL_DETERMINANT(g11(i,j,k), g12(i,j,k), g13(i,j,k), \ + g22(i,j,k), g23(i,j,k), g33(i,j,k)) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 @@ -347,8 +456,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),gaa(i,j,k),gab(i,j,k),& - gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),g12(i,j,k),& + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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 +469,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, & - 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, & + g11(i,j,k), g12(i,j,k), g13(i,j,k), & + g22(i,j,k), g23(i,j,k), g33(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 +483,8 @@ subroutine GRHydro_InitialAtmosphereReset(CCTK_ARGUMENTS) vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0 - 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)) + det = SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k)) if(evolve_temper.ne.0) then ! set the temperature to be relatively low @@ -387,8 +496,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),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), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p(i,j,k),g12_p(i,j,k),& + g13_p(i,j,k),g22_p(i,j,k),g23_p(i,j,k),g33_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 +509,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, & - 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, & + g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), & + g22_p(i,j,k), g23_p(i,j,k), g33_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 +524,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(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 = SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k)) if(evolve_temper.ne.0) then ! set the temperature to be relatively low @@ -427,8 +536,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),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), & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p_p(i,j,k),g12_p_p(i,j,k),& + g13_p_p(i,j,k),g22_p_p(i,j,k),g23_p_p(i,j,k),g33_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 +549,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, & - 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, & + g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), & + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_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 6ff530b..414313b 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -56,7 +56,7 @@ SRCS = Utils.F90 \ GRHydro_TmunuM.F90 \ GRHydro_UpdateMaskM.F90 \ GRHydro_UtilsM.F90 \ - GRHydro_TransformTensorBasis.F90 \ + GRHydro_TransformTensorBasis.c \ GRHydro_Jacobian_state.c \ GRHydro_PPMMReconstruct_drv.F90\ GRHydro_ENOReconstruct_drv.F90\ |