aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl13
-rw-r--r--param.ccl35
-rw-r--r--schedule.ccl101
-rw-r--r--src/GRHydro_Con2Prim.F90754
-rw-r--r--src/GRHydro_Con2PrimM.F90283
-rw-r--r--src/GRHydro_HLLE_EOSG.F90576
-rw-r--r--src/GRHydro_HLLE_EOSGM.F90713
-rw-r--r--src/GRHydro_Marquina.F90243
-rw-r--r--src/GRHydro_Prim2Con.F90178
-rw-r--r--src/GRHydro_Prim2ConM.F90166
-rw-r--r--src/GRHydro_Reconstruct.F9018
-rw-r--r--src/GRHydro_ReconstructM.F906
-rw-r--r--src/GRHydro_ReconstructPoly.F9018
-rw-r--r--src/GRHydro_ReconstructPolyM.F909
-rw-r--r--src/GRHydro_RegisterGZ.cc182
-rw-r--r--src/GRHydro_RegisterGZM.cc138
-rw-r--r--src/GRHydro_RiemannSolve.F90231
-rw-r--r--src/GRHydro_RiemannSolveM.F90353
-rw-r--r--src/GRHydro_RoeSolver.F90243
-rw-r--r--src/GRHydro_Scalars.F906
-rw-r--r--src/GRHydro_Startup.F90383
-rw-r--r--src/make.code.defn5
22 files changed, 17 insertions, 4637 deletions
diff --git a/interface.ccl b/interface.ccl
index 953f989..2b4439c 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -289,19 +289,6 @@ void FUNCTION EOS_Omni_EpsFromPress(CCTK_INT IN eoskey, \
USES FUNCTION EOS_Omni_EpsFromPress
-
-
-
-#######################################################
-### Functions provided by the general EOS interface ###
-#######################################################
-
-CCTK_INT FUNCTION EOS_SetupCall(CCTK_INT IN table_handle)
-USES FUNCTION EOS_SetupCall
-
-CCTK_INT FUNCTION EOS_SetGFs(CCTK_POINTER_TO_CONST IN GH, \
- CCTK_INT IN call_number)
-USES FUNCTION EOS_SetGFs
public:
diff --git a/param.ccl b/param.ccl
index d83e028..ac01a95 100644
--- a/param.ccl
+++ b/param.ccl
@@ -479,41 +479,6 @@ keyword gradient_method "Which method is used to set GRHydro::DiffRho?"
"Rho weighted" :: "Based on the size of rho"
} "First diff"
-#For the new EOS interface
-
-BOOLEAN use_eosgeneral "Should we use the new EOS interface?"
-{
-} "no"
-
-STRING eosgeneral_name "The name of the EOS to be used"
-{
- ".*" :: "Any registered EOS"
-} "Polytrope"
-
-STRING eosgeneral_indep_names "The names of additional variables to be used in the EOS call"
-{
- ".*" :: "Any variables consistent with the EOS"
-} ""
-
-INT eosgeneral_n_indeps "The number of additional variables to be used in the EOS call"
-{
- 0:* :: "Anything non-negative"
-} 0
-
-STRING eosgeneral_indep_gfs "The names of the grid functions holding the additional variables"
-{
- ".*" :: "Should be in the full Thorn::GF format"
-} ""
-
-STRING eosgeneral_indep_plus_gfs "The names of the grid functions holding the additional variables that will be correct at the plus cell boundary"
-{
- ".*" :: "Should be in the full Thorn::GF format"
-} ""
-
-STRING eosgeneral_indep_minus_gfs "The names of the grid functions holding the additional variables that will be correct at the minus cell boundary"
-{
- ".*" :: "Should be in the full Thorn::GF format"
-} ""
# NaN detection level
diff --git a/schedule.ccl b/schedule.ccl
index 4e5486a..284c773 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -170,7 +170,7 @@ STORAGE: conformal_state
### Register startup banner ###
###############################
-schedule GRHydro_Startup AT WRAGH AFTER EOSBase_GeneralRegister AFTER HydroBase_StartUp
+schedule GRHydro_Startup AT WRAGH AFTER HydroBase_StartUp
{
LANG: Fortran
} "Startup banner"
@@ -657,26 +657,7 @@ if (CCTK_Equals(method_type, "RSA FV"))
### This is a wrapper call to the Riemann solvers. ###
######################################################
- if (use_eosgeneral)
- {
- if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
- {
- schedule RiemannSolveGeneralM IN FluxTerms AFTER Reconstruct AS Riemann
- {
- LANG: Fortran
- STORAGE: EOS_temps
- STORAGE: RoeAverage_temps
- } "Solve the local Riemann problems - MHD Version"
- } else {
- schedule RiemannSolveGeneral IN FluxTerms AFTER Reconstruct AS Riemann
- {
- LANG: Fortran
- STORAGE: EOS_temps
- STORAGE: RoeAverage_temps
- } "Solve the local Riemann problems"
- }
- }
- else if (CCTK_Equals(GRHydro_eos_type,"General")) {
+ if (CCTK_Equals(GRHydro_eos_type,"General")) {
if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
{
schedule RiemannSolveM IN FluxTerms AFTER Reconstruct AS Riemann
@@ -842,83 +823,7 @@ schedule GRHydro_RefinementLevel AT INITIAL BEFORE HydroBase_Initial
} "Calculate current refinement level"
-if (use_eosgeneral)
-{
-
- STORAGE:Metric_temps
-
- schedule SetMetricTemps IN MoL_PostStep BEFORE HydroBase_PostStep AFTER ADMBase_SetADMVars
- {
- LANG: Fortran
- SYNC: metric_temps
- } "Set the temporary metric terms"
-
- schedule SetMetricTemps AT POSTRESTRICTINITIAL BEFORE HydroBase_PostStep AFTER ADMBase_SetADMVars
- {
- LANG: Fortran
- SYNC: metric_temps
- } "Set the temporary metric terms"
-
- if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
- {
- if (CCTK_Equals(GRHydro_eos_type,"General"))
- {
- schedule Conservative2PrimitiveGeneralM IN HydroBase_Con2Prim AS Con2Prim
- {
- LANG: Fortran
- STORAGE: Con2Prim_temps
- } "Convert back to primitive variables (general) - MHD version"
-
- schedule Primitive2ConservativeCellsGeneralM IN HydroBase_Prim2ConInitial
- {
- LANG: Fortran
- } "Convert initial data given in primive variables to conserved variables - MHD version"
- }
- else if (CCTK_Equals(GRHydro_eos_type,"Polytype"))
- {
- schedule Con2PrimPolytypeGeneralM IN HydroBase_Con2Prim AS Con2Prim
- {
- LANG: Fortran
- STORAGE: Con2Prim_temps
- } "Convert back to primitive variables (polytype) - MHD version"
-
- schedule Primitive2ConservativeCellsGeneralM IN HydroBase_Prim2ConInitial
- {
- LANG: Fortran
- } "Convert initial data given in primive variables to conserved variables - MHD version"
- }
- } else {
- if (CCTK_Equals(GRHydro_eos_type,"General"))
- {
- schedule Conservative2PrimitiveGeneral IN HydroBase_Con2Prim AS Con2Prim
- {
- LANG: Fortran
- STORAGE: Con2Prim_temps
- } "Convert back to primitive variables (general)"
-
- schedule Primitive2ConservativeCellsGeneral IN HydroBase_Prim2ConInitial
- {
- LANG: Fortran
- } "Convert initial data given in primive variables to conserved variables"
-
- }
- else if (CCTK_Equals(GRHydro_eos_type,"Polytype"))
- {
- schedule Con2PrimPolytypeGeneral IN HydroBase_Con2Prim AS Con2Prim
- {
- LANG: Fortran
- STORAGE: Con2Prim_temps
- } "Convert back to primitive variables (polytype)"
-
- schedule Primitive2ConservativeCellsGeneral IN HydroBase_Prim2ConInitial
- {
- LANG: Fortran
- } "Convert initial data given in primive variables to conserved variables"
-
- }
- }
-}
-else if (CCTK_Equals(GRHydro_eos_type,"General"))
+if (CCTK_Equals(GRHydro_eos_type,"General"))
{
if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index ae83741..4d0f015 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -1724,760 +1724,6 @@ subroutine Con2Prim_ptBoundsTracer(cons_tracer, tracer, rho, one_over_w_lorentz,
end subroutine Con2Prim_ptBoundsTracer
- /*@@
- @routine Conservative2PrimitiveGeneral
- @date Sat Jan 26 01:08:33 2002
- @author Ian Hawke
- @desc
- Wrapper routine that converts from conserved to primitive variables
- at every grid cell centre. Converted to the EOSGeneral format
- @enddesc
- @calls
- @calledby
- @history
- @endhistory
-
-@@*/
-
-subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, itracer, nx, ny, nz, ierr
- character(len=100) warnline
-
- integer :: maxerrloc(3)
- integer :: ii,jj,kk
- CCTK_REAL :: maxerr
-
-
- CCTK_INT :: count
- CCTK_REAL :: s2, df, v2, w, vlowx, vlowy, vlowz
- CCTK_REAL :: temp1, temp2, temp3
- CCTK_REAL :: udens, usx, usy, usz, utau
- CCTK_REAL :: drhodpress, depsdpress
- logical :: loop_condition
-
- CCTK_REAL, dimension(:,:,:), allocatable :: f
- logical, dimension(:,:,:), allocatable :: point_loop_condition
-
- CCTK_INT :: type_bits, atmosphere
- CCTK_INT :: type2_bits
-
- CCTK_REAL :: local_min_tracer
-
- CCTK_INT, dimension(3) :: loc, loc2
-
- allocate(f(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
- point_loop_condition(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
- STAT=ierr)
- if (ierr .ne. 0) then
- call CCTK_WARN(0, "Allocation problems with f")
- end if
-
- call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
- call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
- type2_bits = -1
-
- nx = cctk_lsh(1)
- ny = cctk_lsh(2)
- nz = cctk_lsh(3)
-
- count = 0
- press_old = max(eosgeneral_pmin, press)
-
- if (use_min_tracer .ne. 0) then
- local_min_tracer = min_tracer
- else
- local_min_tracer = 0d0
- end if
-
- loop_condition = .true.
-
-!!$ Set up rho and epsilon
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- if (evolve_tracer .ne. 0) then
- do itracer=1,number_of_tracers
- call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), &
- tracer(i,j,k,itracer), dens(i,j,k))
- if (use_min_tracer .ne. 0) then
- if (tracer(i,j,k,itracer) .le. local_min_tracer) then
- tracer(i,j,k,itracer) = local_min_tracer
- end if
- end if
- enddo
- endif
-
- if(evolve_Y_e.ne.0) then
- Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
- endif
-
- if ( (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
- GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) &
- .or.(tau(i,j,k) .le. 0d0)) then
-
- rho(i,j,k) = GRHydro_rho_min
- dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
- scon(i,j,k,:) = 0.d0
- w_lorentz(i,j,k) = 1.d0
- press(i,j,k) = 0.1d0 * eosgeneral_pmin
-! eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved
- eps(i,j,k) = 1.0e-12
- tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * eps(i,j,k)
-
- SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
- atmosphere_mask(i,j,k) = 1
-
- end if
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
- s2 = usx*usx*GRHydro_uxx(i,j,k) + &
- usy*usy*GRHydro_uyy(i,j,k) + &
- usz*usz*GRHydro_uzz(i,j,k) + &
- 2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
- 2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
- 2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
-!!$ Check that the variables have a chance of being physical
-
- if( (utau + press_old(i,j,k) + udens)**2 - s2 .le. 0.0d0) then
- if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
- !$OMP CRITICAL
- call CCTK_WARN(1,'Con2PrimGeneral: variables unphysical!')
- write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,g15.6)') 'undensitized energy: ',utau
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,g15.6)') 'undensitized density: ',udens
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,g15.6)') 'undensitized s2: ',s2
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,g15.6)') 'pressure guess: ',press_old(i,j,k)
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,i4)') 'refinement level: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
- x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
- call CCTK_WARN(1,warnline)
- call CCTK_WARN(GRHydro_c2p_warnlevel, "Unphysical variables")
- !$OMP END CRITICAL
- exit
- else
- !$OMP CRITICAL
- write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
- x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
- call CCTK_WARN(1,warnline)
- !$OMP END CRITICAL
- exit
- endif
- endif
-
- temp1 = (utau + udens + press_old(i,j,k))
- temp2 = temp1**2 - s2
- temp3 = sqrt(temp2)
-
- rho(i,j,k) = udens * temp3 / temp1
- w_lorentz(i,j,k) = temp1 / temp3
-
-!!$ Eps not previously set?
- eps(i,j,k) = &
- (temp3 - press_old(i,j,k) * w_lorentz(i,j,k) - udens) / &
- udens
-
- end do
- end do
- end do
-
- ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
- where (press .le. eosgeneral_pmin)
- press = eosgeneral_pmin
- rho = GRHydro_rho_min
- dens = sqrt(GRHydro_det) * rho
- scon(:,:,:,1) = 0.d0
- scon(:,:,:,2) = 0.d0
- scon(:,:,:,3) = 0.d0
- w_lorentz = 1d0
-! eps = press / rho
- eps = 1.0e-12
- tau = sqrt(GRHydro_det) * rho * eps
- end where
- f = press_old - press
- press_new = press_old
-
- do while (loop_condition)
-
- count = count + 1
-
- if (count > GRHydro_countmax) then
- if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
- !$OMP CRITICAL
- call CCTK_WARN(1, 'Con2PrimGeneral: error: did not converge in ')
- write(warnline,'(a20,i12,a10)') ' ',&
- GRHydro_countmax,' steps'
- call CCTK_WARN(1,warnline)
- write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- maxerr = maxval(abs(press_new - press_old) &
- / abs(press_new))
- maxerrloc = maxloc(abs(press_new - press_old) &
- / abs(press_new))
- ii = maxerrloc(1)
- jj = maxerrloc(2)
- kk = maxerrloc(3)
- write(warnline,'(1P10E15.6)') maxerr
- call CCTK_WARN(1,warnline)
- write(warnline,'(1P10E15.6)') press_new(ii,jj,kk), press_old(ii,jj,kk)
- call CCTK_WARN(1,warnline)
- write(warnline,'(3i4)') ii,jj,kk
- call CCTK_WARN(1,warnline)
- write(warnline,'(1P10E15.6)') x(ii,jj,kk), &
- y(ii,jj,kk), &
- z(ii,jj,kk)
- call CCTK_WARN(1,warnline)
- call CCTK_WARN(0,"Abort!")
- !$OMP END CRITICAL
- exit
- else
- !$OMP CRITICAL
- write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',&
- GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- !$OMP END CRITICAL
- exit
- endif
- end if
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- !do not compute if in atmosphere or in an excised region
- 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)) &
- then
- press_new(i,j,k) = press(i,j,k)
- press_old(i,j,k) = press(i,j,k)
-! cycle
- endif
-
-
- ! All right, is this point already converged?!?
- loop_condition = abs(press_new(i,j,k)-press_old(i,j,k)) / abs(press_new(i,j,k)) > &
- GRHydro_perc_ptol
-
- loop_condition = loop_condition .and. &
- abs(press_new(i,j,k)-press_old(i,j,k)) > GRHydro_del_ptol
-
- loop_condition = loop_condition .or. (count < GRHydro_countmin) &
- .or. (count < 2) ! we have to at least go once through c2p.
-
- point_loop_condition(i,j,k) = loop_condition
-
- if(.not.loop_condition) cycle
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
- s2 = usx*usx*GRHydro_uxx(i,j,k) + &
- usy*usy*GRHydro_uyy(i,j,k) + &
- usz*usz*GRHydro_uzz(i,j,k) + &
- 2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
- 2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
- 2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
- temp1 = (utau + udens + press_new(i,j,k))
- temp2 = temp1**2 - s2
- temp3 = sqrt(temp2)
-
- drhodpress = udens * s2 / (temp3 * temp1**2)
- depsdpress = press_new(i,j,k) * s2 / (rho(i,j,k) * temp1 * temp2)
-
- df = 1.d0 - eos_dpdeps_temp(i,j,k) * depsdpress - &
- eos_dpdrho_temp(i,j,k) * drhodpress
-
- press_old(i,j,k) = press_new(i,j,k)
-
- press_new(i,j,k) = press_old(i,j,k) - f(i,j,k) / df
- if (press_new(i,j,k) .le. eosgeneral_pmin) then
- press_new(i,j,k) = eosgeneral_pmin
- rho(i,j,k) = GRHydro_rho_min
- dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
- scon(i,j,k,:) = 0.d0
- w_lorentz(i,j,k) = 1d0
- eps(i,j,k) = press(i,j,k) / rho(i,j,k)
- tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * eps(i,j,k)
- udens = rho(i,j,k)
- utau = rho(i,j,k) * eps(i,j,k)
- s2 = 0.d0
- end if
-
- temp1 = (utau + udens + press_new(i,j,k))
- temp2 = temp1**2 - s2
- temp3 = sqrt(temp2)
-
- rho(i,j,k) = udens * temp3 / temp1
- w_lorentz(i,j,k) = temp1 / temp3
- eps(i,j,k) = &
- (temp3 - press_new(i,j,k) * w_lorentz(i,j,k) - udens) / &
- udens
-
- point_loop_condition(i,j,k) = &
- (abs(press_new(i,j,k) - press_old(i,j,k)) / abs(press_new(i,j,k)) > &
- GRHydro_perc_ptol)
- point_loop_condition(i,j,k) = point_loop_condition(i,j,k) .and. &
- (abs(press_new(i,j,k) - press_old(i,j,k)) > &
- GRHydro_del_ptol)
-
- end do
- end do
- end do
-
- ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
- press = max(press, eosgeneral_pmin)
-
- f = press_new - press
-
- loop_condition = any(point_loop_condition)
-
- loop_condition = loop_condition .or. &
- (count < GRHydro_countmin)
-
- end do
-
- deallocate(point_loop_condition)
-
-
- do count = 1, GRHydro_polish
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- !do not compute if in atmosphere or in an excised region
- 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
-
- f(i,j,k) = press_new(i,j,k) - press(i,j,k)
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
- s2 = usx*usx*GRHydro_uxx(i,j,k) + &
- usy*usy*GRHydro_uyy(i,j,k) + &
- usz*usz*GRHydro_uzz(i,j,k) + &
- 2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
- 2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
- 2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
- temp1 = (utau + udens + press_new(i,j,k))
- temp2 = temp1**2 - s2
- temp3 = sqrt(temp2)
-
- drhodpress = udens * s2 / (temp3 * temp1**2)
- depsdpress = press_new(i,j,k) * s2 / (rho(i,j,k) * temp1 * temp2)
-
- df = 1.d0 - eos_dpdeps_temp(i,j,k) * depsdpress - &
- eos_dpdrho_temp(i,j,k) * drhodpress
-
- press_new(i,j,k) = max(press_new(i,j,k) - f(i,j,k) / df, eosgeneral_pmin)
-
- temp1 = (utau + udens + press_new(i,j,k))
- temp2 = temp1**2 - s2
- temp3 = sqrt(temp2)
-
- rho(i,j,k) = udens * temp3 / temp1
- w_lorentz(i,j,k) = temp1 / temp3
- eps(i,j,k) = &
- (temp3 - press_new(i,j,k) * w_lorentz(i,j,k) - udens) / &
- udens
-
- end do
- end do
- end do
-
- ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
- press = max(press, eosgeneral_pmin)
-
- end do
-
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- !do not compute if in atmosphere or in an excised region
- 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
-
-!!$ Note that in the following we enforce eps > 0
-!!$ This makes the later check redundant
-
- if ( (rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) .or. &
- (eps(i,j,k) .lt. 0d0) ) then
-
- rho(i,j,k) = GRHydro_rho_min
- udens = rho(i,j,k)
- utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
- dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
- eps(i,j,k) = utau / udens
- scon(i,j,k,:) = 0.d0
- w_lorentz(i,j,k) = 1.d0
- usx = 0.d0
- usy = 0.d0
- usz = 0.d0
-
- else
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
-
- end if
-
- vlowx = usx / &
- ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
- w_lorentz(i,j,k)**2)
- vlowy = usy / &
- ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
- w_lorentz(i,j,k)**2)
- vlowz = usz / &
- ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
- w_lorentz(i,j,k)**2)
- vel(i,j,k,1) = GRHydro_uxx(i,j,k) * vlowx + &
- GRHydro_uxy(i,j,k) * vlowy + &
- GRHydro_uxz(i,j,k) * vlowz
- vel(i,j,k,2) = GRHydro_uxy(i,j,k) * vlowx + &
- GRHydro_uyy(i,j,k) * vlowy + &
- GRHydro_uyz(i,j,k) * vlowz
- vel(i,j,k,3) = GRHydro_uxz(i,j,k) * vlowx + &
- GRHydro_uyz(i,j,k) * vlowy + &
- GRHydro_uzz(i,j,k) * vlowz
-
-!!$ If eps negative then complain vociferously
-
- if (eps(i,j,k) .le. 0.d0) then
- if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
- !$OMP CRITICAL
- call CCTK_WARN(1,'Con2PrimGeneral: stopping the code.')
- call CCTK_WARN(1, ' specific internal energy just went below 0! ')
- write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- write(warnline,'(a20,3i6)') 'ijk location: ',&
- i,j,k
- call CCTK_WARN(1,warnline)
- write(warnline,'(a20,1P10E15.6)') 'xyz location: ',&
- x(i,j,k),y(i,j,k),z(i,j,k)
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,i5)') 'refinement level: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- write(warnline,'(a20,1P1E15.6)') 'radius: ',r(i,j,k)
- call CCTK_WARN(1,warnline)
- write(warnline,'(a20,1P1E15.6)') 'velocities: ',&
- vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3)
- call CCTK_WARN(1,warnline)
- call CCTK_WARN(GRHydro_c2p_warnlevel, \
- "Specific internal energy negative")
- !$OMP END CRITICAL
- exit
- else
- !$OMP CRITICAL
- write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but \
- I was told to ignore level: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
- x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
- call CCTK_WARN(1,warnline)
- !$OMP END CRITICAL
- exit
- endif
- endif
-
- end do
- end do
- end do
-
- deallocate(f)
-
-end subroutine Conservative2PrimitiveGeneral
-
-
- /*@@
- @routine Conservative2PrimitivePolytypeGeneral
- @date Sat Jan 26 01:08:33 2002
- @author Ian Hawke
- @desc
- Wrapper routine that converts from conserved to primitive variables
- at every grid cell centre. Converted to the EOSGeneral format
- @enddesc
- @calls
- @calledby
- @history
- @endhistory
-
-@@*/
-
-subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, itracer, nx, ny, nz
- character(len=100) warnline
-
- CCTK_INT :: count, ierr
- CCTK_REAL :: s2, f, df, v2, w, vlowx, vlowy, vlowz
- CCTK_REAL :: temp1, temp2, temp3
- CCTK_REAL :: udens, usx, usy, usz
- CCTK_REAL :: enthalpy, denthalpy
- logical :: loop_condition
-
- CCTK_INT :: type_bits
- CCTK_INT :: atmosphere
- CCTK_INT :: type2_bits
-
- CCTK_REAL :: local_min_tracer
-
- call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
- call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
- type2_bits = -1
-
- nx = cctk_lsh(1)
- ny = cctk_lsh(2)
- nz = cctk_lsh(3)
-
- if (use_min_tracer .ne. 0) then
- local_min_tracer = min_tracer
- else
- local_min_tracer = 0d0
- end if
-
-!!$ In what follows press_temp is really rho_temp
-
- count = 0
- press_new = max(GRHydro_rho_min, rho)
-
- loop_condition = .true.
-
- ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-
-!!$ Set up W and epsilon
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- if (evolve_tracer .ne. 0) then
- do itracer=1,number_of_tracers
- call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), &
- tracer(i,j,k,itracer), dens(i,j,k))
-
- if (use_min_tracer .ne. 0) then
- if (tracer(i,j,k,itracer) .le. local_min_tracer) then
- tracer(i,j,k,itracer) = local_min_tracer
- end if
- end if
-
- enddo
- end if
-
- if(evolve_Y_e.ne.0) then
- Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
- endif
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- s2 = usx*usx*GRHydro_uxx(i,j,k) + &
- usy*usy*GRHydro_uyy(i,j,k) + &
- usz*usz*GRHydro_uzz(i,j,k) + &
- 2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
- 2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
- 2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
- enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
- w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) )
-
- end do
- end do
- end do
-
- do while (loop_condition)
-
- count = count + 1
-
- if (count > GRHydro_countmax) then
- if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
- !$OMP CRITICAL
- call CCTK_WARN(1, 'Con2PrimPolytypeGeneral: error: did not converge in ')
- write(warnline,'(a20,g12.7,a10)') ' ',&
- GRHydro_countmax,' steps'
- call CCTK_WARN(1,warnline)
- write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- call CCTK_WARN(GRHydro_c2p_warnlevel, "Did not converge")
- !$OMP END CRITICAL
- exit
- else
- !$OMP CRITICAL
- write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',GRHydro_reflevel
- call CCTK_WARN(1,warnline)
- !$OMP END CRITICAL
- exit
- endif
- end if
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- !do not compute if in atmosphere or in an excised region
- 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
-
- enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
-
- tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * &
- ( press_new(i,j,k) * enthalpy * w_lorentz(i,j,k)**2 - &
- press(i,j,k) ) - dens(i,j,k)
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- s2 = usx*usx*GRHydro_uxx(i,j,k) + &
- usy*usy*GRHydro_uyy(i,j,k) + &
- usz*usz*GRHydro_uzz(i,j,k) + &
- 2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
- 2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
- 2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
- f = press_new(i,j,k) * w_lorentz(i,j,k) - udens
-
- denthalpy = eos_dpdrho_temp(i,j,k) / press_new(i,j,k)
-
- df = w_lorentz(i,j,k) - press_new(i,j,k) * s2 * denthalpy / &
- ( w_lorentz(i,j,k) * (udens**2) * (enthalpy**3) )
-
- rho(i,j,k) = press_new(i,j,k)
- press_new(i,j,k) = max(press_new(i,j,k) - f / df, GRHydro_rho_min)
-
- enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
- w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) )
-
- end do
- end do
- end do
-
- ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-
- loop_condition = (maxval(abs(rho - press_new) / abs(press_new)) > &
- GRHydro_perc_ptol)
- loop_condition = loop_condition .and. &
- (maxval(abs(rho - press_new)) > GRHydro_del_ptol)
- loop_condition = loop_condition .or. &
- (count < GRHydro_countmin)
-
- end do
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- !do not compute if in atmosphere or in an excised region
- 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
-
- udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-
- if (rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
-
- rho(i,j,k) = GRHydro_rho_min
- scon(i,j,k,:) = 0.d0
- w_lorentz(i,j,k) = 1.d0
- usx = 0.d0
- usy = 0.d0
- usz = 0.d0
- s2 = 0.d0
-
- else
-
- usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
- usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
- usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
- s2 = usx*usx*GRHydro_uxx(i,j,k) + &
- usy*usy*GRHydro_uyy(i,j,k) + &
- usz*usz*GRHydro_uzz(i,j,k) + &
- 2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
- 2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
- 2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
- end if
-
- enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / rho(i,j,k)
- w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) )
- tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * &
- ( (rho(i,j,k) * enthalpy) * w_lorentz(i,j,k)**2 - &
- press(i,j,k) ) - dens(i,j,k)
- vlowx = usx / &
- ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
- w_lorentz(i,j,k)**2)
- vlowy = usy / &
- ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
- w_lorentz(i,j,k)**2)
- vlowz = usz / &
- ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
- w_lorentz(i,j,k)**2)
- vel(i,j,k,1) = GRHydro_uxx(i,j,k) * vlowx + &
- GRHydro_uxy(i,j,k) * vlowy + &
- GRHydro_uxz(i,j,k) * vlowz
- vel(i,j,k,2) = GRHydro_uxy(i,j,k) * vlowx + &
- GRHydro_uyy(i,j,k) * vlowy + &
- GRHydro_uyz(i,j,k) * vlowz
- vel(i,j,k,3) = GRHydro_uxz(i,j,k) * vlowx + &
- GRHydro_uyz(i,j,k) * vlowy + &
- GRHydro_uzz(i,j,k) * vlowz
-
-!!$ If eps negative then complain vociferously
-
- if (eps(i,j,k) .lt. 0.0d0) then
- call CCTK_WARN(1, "c2p failed: eps went negative")
- GRHydro_C2P_failed(i,j,k) = 1
- endif
-
- end do
- end do
- end do
-
-end subroutine Con2PrimPolytypeGeneral
-
! subroutines to manage the C2P failure mask
subroutine reset_GRHydro_C2P_failed(CCTK_ARGUMENTS)
diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index 79e24eb..582eec5 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -746,286 +746,3 @@ end subroutine Con2PrimBoundsPolytypeM
!!$ Con2Prim_ptTracer, Con2Prim_BoundsTracer, and Con2Prim_ptBoundsTracer need not be rewritten!
- /*@@
- @routine Conservative2PrimitiveGeneralM
- @date Sep 16, 2010
- @author Scott Noble, Joshua Faber, Bruno Mundim, Ian Hawke
- @desc
- Wrapper routine that converts from conserved to primitive variables
- at every grid cell centre. Converted to the EOSGeneral format
- @enddesc
- @calls
- @calledby
- @history
- @endhistory
-
-@@*/
-
-subroutine Conservative2PrimitiveGeneralM(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
- USE Con2PrimM_fortran_interfaces
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, itracer, nx, ny, nz
- character(len=100) warnline
-
- integer :: maxerrloc(3)
- integer :: ii,jj,kk
- CCTK_REAL :: maxerr
-
-
- CCTK_INT :: count
- CCTK_REAL :: b2,det
- CCTK_REAL :: uxx,uxy,uxz,uyy,uyz,uzz, local_gam, local_pgam
-
- CCTK_INT :: type_bits, atmosphere
- CCTK_INT :: type2_bits,epsnegative
-
- CCTK_REAL :: local_min_tracer
-
- CCTK_INT, dimension(3) :: loc, loc2
-
-!!$ begin EOS Omni vars
- integer :: n,keytemp,anyerr,keyerr(1)
- real*8 :: xpress(1),xtemp(1),xye(1),xdens(1),xeps(1)
- n=1
- keytemp=0
- anyerr=0
- keyerr(1)=0
- xpress(1)=0.0d0
- xtemp(1)=0.0d0
- xye(1)=0.0d0
-!!$ end EOS Omni vars
-
- xdens(1)=1.d0
- xeps(1)=1.d0
- call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
- xdens,xeps,xtemp,xye,xpress,keyerr,anyerr)
- local_gam=xpress(1)+1.d0
-
- call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
- call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
- type2_bits = -1
-
- nx = cctk_lsh(1)
- ny = cctk_lsh(2)
- nz = cctk_lsh(3)
-
- count = 0
- press_old = max(eosgeneral_pmin, press)
-
- if (use_min_tracer .ne. 0) then
- local_min_tracer = min_tracer
- else
- local_min_tracer = 0d0
- end if
-
-!!$ Set up rho and epsilon
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- if (evolve_tracer .ne. 0) then
- do itracer=1,number_of_tracers
- call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), &
- tracer(i,j,k,itracer), dens(i,j,k))
- if (use_min_tracer .ne. 0) then
- if (tracer(i,j,k,itracer) .le. local_min_tracer) then
- tracer(i,j,k,itracer) = local_min_tracer
- end if
- end if
- enddo
- endif
-
- if(evolve_Y_e.ne.0) then
- Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
- endif
-
-
- if ( (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
- GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) &
- .or.(tau(i,j,k) .le. 0d0)) then
-
- b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
- 2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ &
- gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3))
-
-
- rho(i,j,k) = GRHydro_rho_min
- dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
- scon(i,j,k,:) = 0.d0
- w_lorentz(i,j,k) = 1.d0
- press(i,j,k) = 0.1d0 * eosgeneral_pmin
- ! eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved
- eps(i,j,k) = 1.0e-12
- tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * (1.0+eps(i,j,k)+b2/2.0) - dens(i,j,k)
-
- SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
- atmosphere_mask(i,j,k) = 1
-
- end if
-
- det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
- call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
- gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
- gyz(i,j,k),gzz(i,j,k))
-
- call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam,dens(i,j,k), &
- scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k),&
- rho(i,j,k),&
- vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
- gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
- uxx,uxy,uxz,uyy,uyz,uzz,det, &
- Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
- epsnegative,GRHydro_C2P_failed(i,j,k))
-
- end do
- end do
- end do
-
-
-end subroutine Conservative2PrimitiveGeneralM
-
-
- /*@@
- @routine Conservative2PrimitivePolytypeGeneralM
- @date Sep 16, 2010
- @author Scott Noble, Joshua Faber, Bruno Mundim, Ian Hawke
- @desc
- Wrapper routine that converts from conserved to primitive variables
- at every grid cell centre. Converted to the EOSGeneral format
- @enddesc
- @calls
- @calledby
- @history
- @endhistory
-
-@@*/
-
-subroutine Con2PrimPolytypeGeneralM(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, itracer, nx, ny, nz
- character(len=100) warnline
-
- CCTK_REAL :: b2, det
- CCTK_REAL :: uxx,uxy,uxz,uyy,uyz,uzz
-
- CCTK_INT :: type_bits
- CCTK_INT :: atmosphere
- CCTK_INT :: type2_bits,epsnegative
-
- CCTK_REAL :: local_min_tracer, local_pgam, local_k, sc
-
-! begin EOS Omni vars
- integer :: n,keytemp,anyerr,keyerr(1)
- real*8 :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1)
- n=1;keytemp=0;anyerr=0;keyerr(1)=0
- xpress(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0;xeps(1)=0.0d0
-! end EOS Omni vars
-
- call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
- call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
- type2_bits = -1
-
- nx = cctk_lsh(1)
- ny = cctk_lsh(2)
- nz = cctk_lsh(3)
-
- if (use_min_tracer .ne. 0) then
- local_min_tracer = min_tracer
- else
- local_min_tracer = 0d0
- end if
-
-!!$ In what follows press_temp is really rho_temp
-
- do k = 1, nz
- do j = 1, ny
- do i = 1, nx
-
- if (evolve_tracer .ne. 0) then
- do itracer=1,number_of_tracers
- call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), &
- tracer(i,j,k,itracer), dens(i,j,k))
-
- if (use_min_tracer .ne. 0) then
- if (tracer(i,j,k,itracer) .le. local_min_tracer) then
- tracer(i,j,k,itracer) = local_min_tracer
- end if
- end if
- enddo
- end if
-
- if(evolve_Y_e.ne.0) then
- Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
- endif
-
-
- if (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
- GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) then
-
- b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
- 2.0*(gxy(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,2)+gxz(i,j,k)*Bvec(i,j,k,1)*Bvec(i,j,k,3)+ &
- gyz(i,j,k)*Bvec(i,j,k,2)*Bvec(i,j,k,3))
-
-
- rho(i,j,k) = GRHydro_rho_min
- dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
- scon(i,j,k,:) = 0.d0
- w_lorentz(i,j,k) = 1.d0
- press(i,j,k) = 0.1d0 * eosgeneral_pmin
- ! eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved
- eps(i,j,k) = 1.0e-12
- tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * (1.0+eps(i,j,k)+b2/2.0) - dens(i,j,k)
-
- SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
- atmosphere_mask(i,j,k) = 1
-
- end if
-
- det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
- call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,&
- gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
- gyz(i,j,k),gzz(i,j,k))
-
- xrho(1)=1.d0
- call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
- xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
- local_K = xpress(1)
-
- xrho(1)=10.d0
- call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
- xrho,xeps,xtemp,xye,xpress,keyerr,anyerr)
- local_pgam=log(xpress(1)/local_K)/log(xrho(1))
- sc = local_K*dens(i,j,k)
-
- call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,dens(i,j,k), &
- scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc,&
- rho(i,j,k),&
- vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
- gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
- uxx,uxy,uxz,uyy,uyz,uzz,det, &
- Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
- epsnegative,GRHydro_C2P_failed(i,j,k))
-
- end do
- end do
-end do
-
-
-end subroutine Con2PrimPolytypeGeneralM
diff --git a/src/GRHydro_HLLE_EOSG.F90 b/src/GRHydro_HLLE_EOSG.F90
deleted file mode 100644
index 7a3f5f8..0000000
--- a/src/GRHydro_HLLE_EOSG.F90
+++ /dev/null
@@ -1,576 +0,0 @@
- /*@@
- @file GRHydro_HLLE.F90
- @date Sat Jan 26 01:40:14 2002
- @author Ian Hawke, Pedro Montero, Toni Font
- @desc
- The HLLE solver. Called from the wrapper function, so works in
- all directions.
- @enddesc
- @@*/
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-#include "cctk_Functions.h"
-#include "GRHydro_Macros.h"
-#include "SpaceMask.h"
-
- /*@@
- @routine GRHydro_HLLEGeneral
- @date Sat Jan 26 01:41:02 2002
- @author Ian Hawke, Pedro Montero, Toni Font
- @desc
- The HLLE solver. Sufficiently simple that its just one big routine.
- Rewritten for the new EOS interface.
- @enddesc
- @calls
- @calledby
- @history
- Altered from Cactus 3 routines originally written by Toni Font.
- @endhistory
-
-@@*/
-
-subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Eigenproblem
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- integer :: i, j, k, m
- CCTK_REAL, dimension(5) :: cons_p,cons_m,fplus,fminus,lamplus
- CCTK_REAL, dimension(5) :: f1,lamminus
- CCTK_REAL, dimension(5) :: qdiff
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det
- CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
- uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, &
- cs2_p, cs2_m, dpdeps_p, dpdeps_m
-
- CCTK_INT :: type_bits, trivial, not_trivial
-
- integer tadmor
-
- if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
- tadmor = 1
- else
- tadmor = 0
- endif
-
-
- if (flux_direction == 1) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
- &"trivial")
- else if (flux_direction == 2) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
- &"trivial")
- else if (flux_direction == 3) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
- &"trivial")
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
- f1 = 0.d0
- lamminus = 0.d0
- lamplus = 0.d0
- cons_p = 0.d0
- cons_m = 0.d0
- fplus = 0.d0
- fminus = 0.d0
- qdiff = 0.d0
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(6) = pressplus(i,j,k)
- cs2_p = eos_cs2_p(i,j,k)
- dpdeps_p = eos_dpdeps_p(i,j,k)
-
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
- cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
- dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset)
-
-!!$ Calculate various metric terms here.
-!!$ Note also need the average of the lapse at the
-!!$ left and right points.
-
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
-!!$ If the Riemann problem is trivial, just calculate the fluxes from the
-!!$ left state and skip to the next cell
-
- if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then
-
- if (flux_direction == 1) then
- call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
- f1(1),f1(2),f1(3),f1(4),f1(5),&
- prim_m(2),prim_m(6), &
- avg_det,avg_alp,avg_beta)
- else if (flux_direction == 2) then
- call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
- f1(1),f1(3),f1(4),f1(2),f1(5),&
- prim_m(3),prim_m(6), &
- avg_det,avg_alp,avg_beta)
- else if (flux_direction == 3) then
- call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
- f1(1),f1(4),f1(2),f1(3),f1(5), &
- prim_m(4),prim_m(6), &
- avg_det,avg_alp,avg_beta)
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- else !!! The end of this branch is right at the bottom of the routine
-
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh,gxyh,gxzh, &
- gyyh,gyzh, gzzh)
-
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
-!!$ Calculate the jumps in the conserved variables
-
- qdiff(1) = cons_m(1) - cons_p(1)
- qdiff(2) = cons_m(2) - cons_p(2)
- qdiff(3) = cons_m(3) - cons_p(3)
- qdiff(4) = cons_m(4) - cons_p(4)
- qdiff(5) = cons_m(5) - cons_p(5)
-
-!!$ Eigenvalues and fluxes either side of the cell interface
-
- if (flux_direction == 1) then
- call eigenvalues_general(&
- prim_m(2),prim_m(3),prim_m(4),cs2_m, &
- lamminus,&
- gxxh,gxyh,gxzh,&
- gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_general(&
- prim_p(2),prim_p(3),prim_p(4),cs2_p, &
- lamplus,&
- gxxh,gxyh,gxzh,&
- gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
- fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),&
- prim_p(2),prim_p(6), &
- avg_det,avg_alp,avg_beta)
- call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
- fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),&
- prim_m(2),prim_m(6), &
- avg_det,avg_alp,avg_beta)
- else if (flux_direction == 2) then
- call eigenvalues_general(&
- prim_m(3),prim_m(4),prim_m(2),cs2_m, &
- lamminus,&
- gyyh,gyzh,gxyh,&
- gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_general(&
- prim_p(3),prim_p(4),prim_p(2),cs2_p, &
- lamplus,&
- gyyh,gyzh,gxyh,&
- gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
- fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),&
- prim_p(3),prim_p(6), &
- avg_det,avg_alp,avg_beta)
- call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
- fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),&
- prim_m(3),prim_m(6), &
- avg_det,avg_alp,avg_beta)
- else if (flux_direction == 3) then
- call eigenvalues_general(&
- prim_m(4),prim_m(2),prim_m(3),cs2_m, &
- lamminus,&
- gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_general(&
- prim_p(4),prim_p(2),prim_p(3),cs2_p, &
- lamplus,&
- gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
- fplus(1),fplus(4),fplus(2),fplus(3),fplus(5), &
- prim_p(4),prim_p(6), &
- avg_det,avg_alp,avg_beta)
- call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
- fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), &
- prim_m(4),prim_m(6), &
- avg_det,avg_alp,avg_beta)
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- if(tadmor.eq.0) then
-
-!!$ Find minimum and maximum wavespeeds
-
- charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
- charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
- charpm = charmax - charmin
-
-!!$ Calculate flux by standard formula
-
- do m = 1,5
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
-
- end do
-
- else
- ! Tadmor's semi-discrete scheme: JcP 160, 241 (2000)
-
- charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), &
- abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), &
- abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5)))
-
- do m = 1,5
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
- qdiff(m)
-
- end do
-
- end if
-
-
- end if !!! The end of the SpaceMask check for a trivial RP.
-
- densflux(i, j, k) = f1(1)
- sxflux(i, j, k) = f1(2)
- syflux(i, j, k) = f1(3)
- szflux(i, j, k) = f1(4)
- tauflux(i, j, k) = f1(5)
-
- if(evolve_Y_e.ne.0) then
- if (densflux(i, j, k) > 0.d0) then
- Y_e_con_flux(i, j, k) = &
- Y_e_plus(i, j, k) * &
- densflux(i, j, k)
- else
- Y_e_con_flux(i, j, k) = &
- Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * &
- densflux(i, j, k)
- endif
- endif
-
- end do
- end do
- end do
-
-
-end subroutine GRHydro_HLLEGeneral
-
-
-subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Eigenproblem
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- integer :: i, j, k, m
- CCTK_REAL, dimension(number_of_tracers) :: cons_p,cons_m,fplus,fminus,f1
- CCTK_REAL, dimension(5) :: lamminus,lamplus
- CCTK_REAL, dimension(number_of_tracers) :: qdiff
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det
- CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
- uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, &
- cs2_p, cs2_m, dpdeps_p, dpdeps_m
-
- integer tadmor
-
- if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
- tadmor = 1
- else
- tadmor = 0
- endif
-
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
- f1 = 0.d0
- lamminus = 0.d0
- lamplus = 0.d0
- cons_p = 0.d0
- cons_m = 0.d0
- fplus = 0.d0
- fminus = 0.d0
- qdiff = 0.d0
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(:) = cons_tracerplus(i,j,k,:)
- cons_m(:) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(6) = pressplus(i,j,k)
- cs2_p = eos_cs2_p(i,j,k)
- dpdeps_p = eos_dpdeps_p(i,j,k)
-
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
- cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
- dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset)
-
-!!$ Calculate various metric terms here.
-!!$ Note also need the average of the lapse at the
-!!$ left and right points.
-
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
-!!$ If the Riemann problem is trivial, just calculate the fluxes from the
-!!$ left state and skip to the next cell
-
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, &
- gyyh, gyzh, gzzh)
-
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
-!!$ Eigenvalues and fluxes either side of the cell interface
-
- if (flux_direction == 1) then
- call eigenvalues_general(&
- prim_m(2),prim_m(3),prim_m(4),cs2_m, &
- lamminus,&
- gxxh,gxyh,gxzh,&
- gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_general(&
- prim_p(2),prim_p(3),prim_p(4),cs2_p, &
- lamplus,&
- gxxh,gxyh,gxzh,&
- gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * &
- cons_tracerplus(i,j,k,:)
- fminus(:) = (velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
- cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
- else if (flux_direction == 2) then
- call eigenvalues_general(&
- prim_m(3),prim_m(4),prim_m(2),cs2_m, &
- lamminus,&
- gyyh,gyzh,gxyh,&
- gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_general(&
- prim_p(3),prim_p(4),prim_p(2),cs2_p, &
- lamplus,&
- gyyh,gyzh,gxyh,&
- gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * &
- cons_tracerplus(i,j,k,:)
- fminus(:) = (velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
- cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
- else if (flux_direction == 3) then
- call eigenvalues_general(&
- prim_m(4),prim_m(2),prim_m(3),cs2_m, &
- lamminus,&
- gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_general(&
- prim_p(4),prim_p(2),prim_p(3),cs2_p, &
- lamplus,&
- gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * &
- cons_tracerplus(i,j,k,:)
- fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
- cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- if(tadmor.eq.0) then
-
-!!$ Find minimum and maximum wavespeeds
-
- charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
- charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
-
- charpm = charmax - charmin
-
-!!$ Calculate flux by standard formula
-
- do m = 1,number_of_tracers
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
- end do
-
- else
- ! Tadmor's semi-descrite scheme: JcP 160, 241 (2000)
-
- charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), &
- abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), &
- abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5)))
-
- do m = 1,number_of_tracers
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
- qdiff(m)
-
- end do
-
-
- end if
-
-
- cons_tracerflux(i,j,k,:) = f1(:)
-
- end do
- end do
- end do
-
-end subroutine GRHydro_HLLE_TracerGeneral
-
-
-
-
diff --git a/src/GRHydro_HLLE_EOSGM.F90 b/src/GRHydro_HLLE_EOSGM.F90
deleted file mode 100644
index c89d047..0000000
--- a/src/GRHydro_HLLE_EOSGM.F90
+++ /dev/null
@@ -1,713 +0,0 @@
- /*@@
- @file GRHydro_HLLEM.F90
- @date Aug 30, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke, Pedro Montero, Toni Font
- @desc
- The HLLE solver. Called from the wrapper function, so works in
- all directions.
- @enddesc
- @@*/
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-#include "cctk_Functions.h"
-#include "GRHydro_Macros.h"
-#include "SpaceMask.h"
-
- /*@@
- @routine GRHydro_HLLEGeneralM
- @date Aug 30, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke, Pedro Montero, Toni Font
- @desc
- The HLLE solver. Sufficiently simple that its just one big routine.
- Rewritten for the new EOS interface.
- @enddesc
- @calls
- @calledby
- @history
- Altered from Cactus 3 routines originally written by Toni Font.
- @endhistory
-
-@@*/
-
-subroutine GRHydro_HLLEGeneralM(CCTK_ARGUMENTS)
- USE GRHydro_EigenproblemM
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- integer :: i, j, k, m
- CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL, dimension(5) :: lamminus,lamplus
- CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det
- CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
- uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, &
- cs2_p, cs2_m, dpdeps_p, dpdeps_m
- CCTK_REAL :: rhoenth_p, rhoenth_m, avg_betax, avg_betay, avg_betaz
- CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm
- CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p
- CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm
- CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm
-
- CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm
-
- CCTK_INT :: type_bits, trivial, not_trivial
-
- integer tadmor
-
- if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
- tadmor = 1
- else
- tadmor = 0
- endif
-
-
- if (flux_direction == 1) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
- &"trivial")
- else if (flux_direction == 2) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
- &"trivial")
- else if (flux_direction == 3) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
- &"trivial")
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
- f1 = 0.d0
- lamminus = 0.d0
- lamplus = 0.d0
- cons_p = 0.d0
- cons_m = 0.d0
- fplus = 0.d0
- fminus = 0.d0
- qdiff = 0.d0
-
- if(clean_divergence.ne.0) then
- psidcp = 0.d0
- psidcm = 0.d0
- endif
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
- cons_p(6) = Bvecxplus(i,j,k)
- cons_p(7) = Bvecyplus(i,j,k)
- cons_p(8) = Bveczplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(6) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(7) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(8) = Bveczminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(6) = pressplus(i,j,k)
- cs2_p = eos_cs2_p(i,j,k)
- dpdeps_p = eos_dpdeps_p(i,j,k)
- rhoenth_p = prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)
-
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
- cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
- dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset)
- rhoenth_m = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)
-
- if(clean_divergence.ne.0) then
- psidcp = psidcplus(i,j,k)
- psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset)
- endif
-
-!!$ Calculate various metric terms here.
-!!$ Note also need the average of the lapse at the
-!!$ left and right points.
-!!$
-!!$ In MHD, we need all three shift components regardless of the flux direction
-
- if (shift_state .ne. 0) then
- avg_betax = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- avg_betay = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- avg_betaz = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- if (flux_direction == 1) then
- avg_beta=avg_betax
- else if (flux_direction == 2) then
- avg_beta=avg_betay
- else if (flux_direction == 3) then
- avg_beta=avg_betaz
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- avg_betax = 0.d0
- avg_betay = 0.d0
- avg_betaz = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
- vxtp = prim_p(2)-avg_betax/avg_alp
- vytp = prim_p(3)-avg_betay/avg_alp
- vztp = prim_p(4)-avg_betaz/avg_alp
- vxtm = prim_m(2)-avg_betax/avg_alp
- vytm = prim_m(3)-avg_betay/avg_alp
- vztm = prim_m(4)-avg_betaz/avg_alp
-
- call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- prim_p(2),prim_p(3),prim_p(4),cons_p(6),cons_p(7),cons_p(8), &
- velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, &
- bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp)
- call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), &
- velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, &
- bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm)
-
- ab0p = wp*bdotvp
- ab0m = wm*bdotvm
-
- vA2p = b2p/(rhoenth_p+b2p)
- vA2m = b2m/(rhoenth_m+b2m)
-
-!!$ p^* = p+b^2/2 in Anton et al.
- pressstarp = prim_p(6)+0.5d0*b2p
- pressstarm = prim_m(6)+0.5d0*b2m
-
-
-!!$ If the Riemann problem is trivial, just calculate the fluxes from the
-!!$ left state and skip to the next cell
-
- if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then
-
-!!$ we now pass in the B-field as conserved and flux, the vtilde's instead of v's,
-!!$ pressstar instead of P, b_i, alp b^0, w, metric determinant,
-!!$ alp, and beta in the flux dir
-
- if (flux_direction == 1) then
- call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
- cons_m(6),cons_m(7),cons_m(8),&
- f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),&
- vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
- if(clean_divergence.ne.0) then
- psidcf = cons_p(6)
- endif
-
- else if (flux_direction == 2) then
- call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
- cons_m(7),cons_m(8),cons_m(6),&
- f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),&
- vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
- if(clean_divergence.ne.0) then
- psidcf = cons_p(7)
- endif
-
- else if (flux_direction == 3) then
- call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
- cons_m(8),cons_m(6),cons_m(7),&
- f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), &
- vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
- if(clean_divergence.ne.0) then
- psidcf = cons_p(8)
- endif
-
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- else !!! The end of this branch is right at the bottom of the routine
-
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, &
- gyyh, gyzh, gzzh)
-
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
-!!$ Calculate the jumps in the conserved variables
-
- qdiff(1) = cons_m(1) - cons_p(1)
- qdiff(2) = cons_m(2) - cons_p(2)
- qdiff(3) = cons_m(3) - cons_p(3)
- qdiff(4) = cons_m(4) - cons_p(4)
- qdiff(5) = cons_m(5) - cons_p(5)
- qdiff(6) = cons_m(6) - cons_p(6)
- qdiff(7) = cons_m(7) - cons_p(7)
- qdiff(8) = cons_m(8) - cons_p(8)
-
- if (clean_divergence.ne.0) then
- psidcdiff = psidcm - psidcp
- endif
-
-!!$ Eigenvalues and fluxes either side of the cell interface
-
- if (flux_direction == 1) then
- call eigenvalues_generalM(&
- prim_m(2),prim_m(3),prim_m(4),cs2_m,vA2m, &
- lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_generalM(&
- prim_p(2),prim_p(3),prim_p(4),cs2_p,vA2p, &
- lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
- cons_p(6),cons_p(7),cons_p(8),&
- fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),&
- vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
- avg_det,avg_alp,avg_beta)
- call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
- cons_m(6),cons_m(7),cons_m(8),&
- fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),&
- fminus(6),fminus(7),fminus(8),&
- vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
-
- if(clean_divergence.ne.0) then
- psidcfp = cons_p(6)
- psidcfm = cons_m(6)
- endif
-
- else if (flux_direction == 2) then
- call eigenvalues_generalM(&
- prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, &
- lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_generalM(&
- prim_p(3),prim_p(4),prim_p(2),cs2_p,vA2p, &
- lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
- cons_p(7),cons_p(8),cons_p(6),&
- fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),&
- vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
- avg_det,avg_alp,avg_beta)
- call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
- cons_m(7),cons_m(8),cons_m(6),&
- fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),&
- fminus(7),fminus(8),fminus(6),&
- vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
- if(clean_divergence.ne.0) then
- psidcfp = cons_p(7)
- psidcfm = cons_m(7)
- endif
-
- else if (flux_direction == 3) then
- call eigenvalues_generalM(&
- prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, &
- lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_generalM(&
- prim_p(4),prim_p(2),prim_p(3),cs2_p,vA2p, &
- lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
- cons_p(8),cons_p(6),cons_p(7),&
- fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), &
- vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
- avg_det,avg_alp,avg_beta)
- call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
- cons_m(8),cons_m(6),cons_m(7),&
- fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), &
- fminus(8),fminus(6),fminus(7), &
- vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
- if(clean_divergence.ne.0) then
- psidcfp = cons_p(8)
- psidcfm = cons_m(8)
- endif
-
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- if(tadmor.eq.0) then
-
-!!$ Find minimum and maximum wavespeeds
-
- charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
- charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
- charpm = charmax - charmin
-
-!!$ Calculate flux by standard formula
-
- do m = 1,8
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
-
- end do
-
- if(clean_divergence.ne.0) then
- psidcdiff = psidcm - psidcp
- psidcf = (charmax * psidcfp - charmin * psidcfm + &
- charmax * charmin * psidcdiff) / charpm
- endif
-
- else
- ! Tadmor's semi-discrete scheme: JcP 160, 241 (2000)
-
- charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), &
- abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), &
- abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5)))
-
- do m = 1,8
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
- qdiff(m)
-
- end do
-
- if(clean_divergence.ne.0) then
- psidcdiff = psidcm - psidcp
- psidcf = 0.5d0 * (psidcfp + psidcfm) - 0.5d0*charmax* &
- psidcdiff
- endif
-
- end if
-
-
-
- end if !!! The end of the SpaceMask check for a trivial RP.
-
- densflux(i, j, k) = f1(1)
- sxflux(i, j, k) = f1(2)
- syflux(i, j, k) = f1(3)
- szflux(i, j, k) = f1(4)
- tauflux(i, j, k) = f1(5)
- Bvecxflux(i, j, k) = f1(6)
- Bvecyflux(i, j, k) = f1(7)
- Bveczflux(i, j, k) = f1(8)
- if(clean_divergence.ne.0) then
- psidcflux(i,j,k) = psidcf
- endif
-
-
- if(evolve_Y_e.ne.0) then
- if (densflux(i, j, k) > 0.d0) then
- Y_e_con_flux(i, j, k) = &
- Y_e_plus(i, j, k) * &
- densflux(i, j, k)
- else
- Y_e_con_flux(i, j, k) = &
- Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * &
- densflux(i, j, k)
- endif
- endif
-
- end do
- end do
- end do
-
-end subroutine GRHydro_HLLEGeneralM
-
-
-subroutine GRHydro_HLLE_TracerGeneralM(CCTK_ARGUMENTS)
-
- USE GRHydro_EigenproblemM
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- integer :: i, j, k, m
- CCTK_REAL, dimension(number_of_tracers) :: cons_p,cons_m,fplus,fminus,f1
- CCTK_REAL, dimension(5) :: lamminus,lamplus
- CCTK_REAL, dimension(number_of_tracers) :: qdiff
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL, dimension(3) :: mag_p, mag_m
- CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det
- CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
- uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, &
- cs2_p, cs2_m, dpdeps_p, dpdeps_m
- CCTK_REAL :: b2p,b2m,vA2m,vA2p
-
- integer tadmor
-
- if(CCTK_EQUALS(HLLE_type,"Tadmor")) then
- tadmor = 1
- else
- tadmor = 0
- endif
-
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
- f1 = 0.d0
- lamminus = 0.d0
- lamplus = 0.d0
- cons_p = 0.d0
- cons_m = 0.d0
- mag_p = 0.d0
- mag_m = 0.d0
- fplus = 0.d0
- fminus = 0.d0
- qdiff = 0.d0
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(:) = cons_tracerplus(i,j,k,:)
- cons_m(:) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
-
- mag_p(1) = Bvecxplus(i,j,k)
- mag_p(2) = Bvecyplus(i,j,k)
- mag_p(3) = Bveczplus(i,j,k)
-
- mag_m(1) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset)
- mag_m(2) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset)
- mag_m(3) = Bveczminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(6) = pressplus(i,j,k)
- cs2_p = eos_cs2_p(i,j,k)
- dpdeps_p = eos_dpdeps_p(i,j,k)
-
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
- cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
- dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset)
-
-!!$ Calculate various metric terms here.
-!!$ Note also need the average of the lapse at the
-!!$ left and right points.
-
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
-!!$ If the Riemann problem is trivial, just calculate the fluxes from the
-!!$ left state and skip to the next cell
-
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, &
- gyyh, gyzh, gzzh)
-
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
-!!$ b^2 = (1-v^2)B^2+(B dot v)^2
- b2p=(1.d0-DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4)))* &
- DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_p(1),mag_p(2),mag_p(3)) + &
- (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4),mag_p(1),mag_p(2),mag_p(3)))**2
- b2m=(1.d0-DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4)))* &
- DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3)) + &
- (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3)))**2
-
- vA2p = b2p/(prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)+b2p)
- vA2m = b2m/(prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)+b2m)
-
-!!$ Eigenvalues and fluxes either side of the cell interface
-
- if (flux_direction == 1) then
- call eigenvalues_generalM(&
- prim_m(2),prim_m(3),prim_m(4),cs2_m,vA2m, &
- lamminus,&
- gxxh,gxyh,gxzh,&
- gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_generalM(&
- prim_p(2),prim_p(3),prim_p(4),cs2_p,vA2p, &
- lamplus,&
- gxxh,gxyh,gxzh,&
- gyyh,gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * &
- cons_tracerplus(i,j,k,:)
- fminus(:) = (velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
- cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
- else if (flux_direction == 2) then
- call eigenvalues_generalM(&
- prim_m(3),prim_m(4),prim_m(2),cs2_m,vA2m, &
- lamminus,&
- gyyh,gyzh,gxyh,&
- gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_generalM(&
- prim_p(3),prim_p(4),prim_p(2),cs2_p,vA2p, &
- lamplus,&
- gyyh,gyzh,gxyh,&
- gzzh,gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * &
- cons_tracerplus(i,j,k,:)
- fminus(:) = (velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
- cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
- else if (flux_direction == 3) then
- call eigenvalues_generalM(&
- prim_m(4),prim_m(2),prim_m(3),cs2_m,vA2m, &
- lamminus,&
- gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues_generalM(&
- prim_p(4),prim_p(2),prim_p(3),cs2_p,vA2p, &
- lamplus,&
- gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * &
- cons_tracerplus(i,j,k,:)
- fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * &
- cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:)
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- if(tadmor.eq.0) then
-
-!!$ Find minimum and maximum wavespeeds
-
- charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
- charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), &
- lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),&
- lamminus(4),lamminus(5))
-
-
- charpm = charmax - charmin
-
-!!$ Calculate flux by standard formula
-
- do m = 1,number_of_tracers
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = (charmax * fplus(m) - charmin * fminus(m) + &
- charmax * charmin * qdiff(m)) / charpm
- end do
-
- else
- ! Tadmor's semi-descrite scheme: JcP 160, 241 (2000)
-
- charmax = max(abs(lamplus(1)), abs(lamplus(2)), abs(lamplus(3)), &
- abs(lamplus(4)),abs(lamplus(5)),abs(lamminus(1)),abs(lamminus(2)), &
- abs(lamminus(3)),abs(lamminus(4)),abs(lamminus(5)))
-
- do m = 1,number_of_tracers
-
- qdiff(m) = cons_m(m) - cons_p(m)
-
- f1(m) = 0.5d0 * (fplus(m) + fminus(m)) - 0.5d0*charmax* &
- qdiff(m)
-
- end do
-
-
- end if
-
-
- cons_tracerflux(i,j,k,:) = f1(:)
-
- end do
- end do
- end do
-
-end subroutine GRHydro_HLLE_TracerGeneralM
-
diff --git a/src/GRHydro_Marquina.F90 b/src/GRHydro_Marquina.F90
index aab367e..8edce62 100644
--- a/src/GRHydro_Marquina.F90
+++ b/src/GRHydro_Marquina.F90
@@ -428,249 +428,6 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS)
return
end subroutine GRHydro_Marquina
- /*@@
- @routine GRHydro_MarquinaGeneral
- @date Wed Feb 13 11:03:32 2002
- @author Pedro Montero, Toni Font
- @desc
- Routine to obtain the Marquina Fluxes
- @enddesc
- @calls
- @calledby
- @history
- Based on routines by Toni Font
- @endhistory
-
-@@*/
-
-
-subroutine GRHydro_MarquinaGeneral(CCTK_ARGUMENTS)
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
- CCTK_REAL, dimension(5) :: cons_p,cons_m,&
- fplus,fminus,marquinaflux
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh, &
- cs2_p,cs2_m,dpdeps_p,dpdeps_m, &
- usendh, psi4h
- integer :: m
- integer :: i,j,k
-
- CCTK_INT :: type_bits, trivial, not_trivial, ierr
-
- if (flux_direction == 1) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
- &"trivial")
- else if (flux_direction == 2) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
- &"trivial")
- else if (flux_direction == 3) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
- &"trivial")
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(6) = pressplus(i,j,k)
- cs2_p = eos_cs2_p(i,j,k)
- dpdeps_p = eos_dpdeps_p(i,j,k)
-
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
- cs2_m = eos_cs2_m(i+xoffset,j+yoffset,k+zoffset)
- dpdeps_m = eos_dpdeps_m(i+xoffset,j+yoffset,k+zoffset)
-
- marquinaflux = 0.d0
-
-!!$ Set metric terms at interface
-
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- !$OMP CRITICAL
- call CCTK_WARN(0, "Flux direction not x,y,z")
- !$OMP END CRITICAL
- end if
- else
- avg_beta = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + &
- gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + &
- gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + &
- gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + &
- gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + &
- gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + &
- gzz(i,j,k))
-
-!!$ routine to calculate the determinant of the metric
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
-!!$ If the Riemann problem is trivial the flux is already correct
-
- if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then
-
- else !!! The end of this branch is right at the bottom of the routine
-
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, gyyh, gyzh, gzzh)
-
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- !$OMP CRITICAL
- call CCTK_WARN(0, "Flux direction not x,y,z")
- !$OMP END CRITICAL
- end if
-
-!!$eigenvalues and right eigenvectors
-
- if (flux_direction == 1) then
-
- call eigenproblem_marquina_general(&
- prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5), &
- prim_m(6),cs2_m,dpdeps_m, &
- prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5), &
- prim_p(6),cs2_p,dpdeps_p, &
- gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- usendh,avg_det,avg_alp,avg_beta, &
- cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5), &
- cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5), &
- marquinaflux(1),marquinaflux(2),marquinaflux(3), &
- marquinaflux(4),marquinaflux(5))
-
- else if (flux_direction == 2) then
-
- call eigenproblem_marquina_general(&
- prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5), &
- prim_m(6),cs2_m,dpdeps_m, &
- prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5), &
- prim_p(6),cs2_p,dpdeps_p, &
- gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
- usendh,avg_det,avg_alp,avg_beta,&
- cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5), &
- cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5), &
- marquinaflux(1),marquinaflux(3),marquinaflux(4), &
- marquinaflux(2),marquinaflux(5))
-
- else if (flux_direction == 3) then
-
- call eigenproblem_marquina_general(&
- prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5), &
- prim_m(6),cs2_m,dpdeps_m, &
- prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5), &
- prim_p(6),cs2_p,dpdeps_p, &
- gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
- usendh,avg_det,avg_alp,avg_beta, &
- cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5), &
- cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5), &
- marquinaflux(1),marquinaflux(4),marquinaflux(2), &
- marquinaflux(3),marquinaflux(5))
-
- else
- !$OMP CRITICAL
- call CCTK_WARN(0, "Flux direction not x,y,z")
- !$OMP END CRITICAL
- end if
-
-!!$ Marquina flux
-
- densflux(i,j,k) = densflux(i,j,k) - 0.5d0 * marquinaflux(1)
- sxflux(i,j,k) = sxflux(i,j,k) - 0.5d0 * marquinaflux(2)
- syflux(i,j,k) = syflux(i,j,k) - 0.5d0 * marquinaflux(3)
- szflux(i,j,k) = szflux(i,j,k) - 0.5d0 * marquinaflux(4)
- tauflux(i,j,k) = tauflux(i,j,k) - 0.5d0 * marquinaflux(5)
-
- end if !!! The end of the SpaceMask check for a trivial RP.
-
- enddo
- enddo
- enddo
-
- if (evolve_tracer .ne. 0) then
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
- if (densflux(i, j, k) > 0.d0) then
-
- cons_tracerflux(i, j, k,:) = &
- tracerplus(i, j, k,:) * &
- densflux(i, j, k)
-
- else
-
- cons_tracerflux(i, j, k,:) = &
- tracerminus(i + xoffset, j + yoffset, k + zoffset,:) * &
- densflux(i, j, k)
-
- end if
-
- end do
- end do
- end do
-
- end if
-
-end subroutine GRHydro_MarquinaGeneral
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90
index 651e228..17548a0 100644
--- a/src/GRHydro_Prim2Con.F90
+++ b/src/GRHydro_Prim2Con.F90
@@ -685,181 +685,3 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS)
end subroutine Prim2ConservativeTracer
- /*@@
- @routine primitive2conservativegeneral
- @date Thu Jan 11 11:03:32 2002
- @author Pedro Montero, Ian Hawke
- @desc
- Converts primitive to conserved variables for the boundary extended data.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
-
-@@*/
-
-subroutine primitive2conservativegeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, ierr
- CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl, &
- gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr, &
- vlowx, vlowy, vlowz
-
- ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallPlus)
- ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallMinus)
-
- !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
- !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr)
- do k = GRHydro_stencil , cctk_lsh(3) - GRHydro_stencil + 1
- do j = GRHydro_stencil , cctk_lsh(2) - GRHydro_stencil + 1
- do i = GRHydro_stencil , cctk_lsh(1) - GRHydro_stencil + 1
-
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
-
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
-
- call prim2conpolytype(-1, gxxl,gxyl,gxzl,&
- gyyl,gyzl,gzzl, &
- 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(-1, gxxr,gxyr,gxzr,&
- gyyr,gyzr,gzzr, &
- 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),&
- velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
- w_lorentzplus(i,j,k))
- end do
- end do
- end do
- !$OMP END PARALLEL DO
-
- ! Note on Y_e: We use Y_e_plus and Y_e_minus directly
- ! in the Riemann solver. That's why it is not necessary
- ! to do a prim2con for Y_e
-
-end subroutine primitive2conservativegeneral
-
- /*@@
- @routine Primitive2ConservativeCellsGeneral
- @date Thu Jan 11 11:03:32 2002
- @author Pedro Montero, Ian Hawke
- @desc
- Converts primitive to conserved variables for the boundary extended data.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
-
-@@*/
-
-subroutine Primitive2ConservativeCellsGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, ierr
- CCTK_REAL :: det,gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt, &
- vlowx, vlowy, vlowz
-
- ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConCellsCall)
-
- do k = 1,cctk_lsh(3)
- do j = 1,cctk_lsh(2)
- do i = 1,cctk_lsh(1)
-
- gxxpt = gxx(i,j,k)
- gxypt = gxy(i,j,k)
- gxzpt = gxz(i,j,k)
- gyypt = gyy(i,j,k)
- gyzpt = gyz(i,j,k)
- gzzpt = gzz(i,j,k)
-
- det = SPATIAL_DETERMINANT(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt)
-
- w_lorentz(i,j,k) = 1.d0 / sqrt(1.d0 - &
- (gxxpt * velx(i,j,k)**2 + &
- gyypt * vely(i,j,k)**2 + &
- gzzpt * velz(i,j,k)**2 + &
- 2.d0 * gxypt * velx(i,j,k) * vely(i,j,k) + &
- 2.d0 * gxzpt * velx(i,j,k) * velz(i,j,k) + &
- 2.d0 * gyzpt * vely(i,j,k) * velz(i,j,k) ))
- vlowx = gxxpt * velx(i,j,k) + &
- gxypt * vely(i,j,k) + &
- gxzpt * velz(i,j,k)
- vlowy = gxypt * velx(i,j,k) + &
- gyypt * vely(i,j,k) + &
- gyzpt * velz(i,j,k)
- vlowz = gxzpt * velx(i,j,k) + &
- gyzpt * vely(i,j,k) + &
- gzzpt * velz(i,j,k)
- dens(i,j,k) = sqrt(det) * rho(i,j,k) * &
- w_lorentz(i,j,k)
- sx(i,j,k) = sqrt(det) * &
- ( rho(i,j,k) * &
- (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
- w_lorentz(i,j,k)**2 * vlowx
- sy(i,j,k) = sqrt(det) * &
- ( rho(i,j,k) * &
- (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
- w_lorentz(i,j,k)**2 * vlowy
- sz(i,j,k) = sqrt(det) * &
- ( rho(i,j,k) * &
- (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
- w_lorentz(i,j,k)**2 * vlowz
- tau(i,j,k) = sqrt(det) * &
- ( ( rho(i,j,k) * &
- (1.d0 + eps(i,j,k)) + press(i,j,k) ) * &
- w_lorentz(i,j,k)**2 - press(i,j,k) ) - &
- dens(i,j,k)
-
- end do
- end do
- end do
-
- if(evolve_Y_e.ne.0) then
- !$OMP PARALLEL DO PRIVATE(i, j)
- 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
- Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k)
- enddo
- enddo
- enddo
- !$OMP END PARALLEL DO
- endif
-
-
-end subroutine Primitive2ConservativeCellsGeneral
diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90
index 90c4644..e50c426 100644
--- a/src/GRHydro_Prim2ConM.F90
+++ b/src/GRHydro_Prim2ConM.F90
@@ -467,169 +467,3 @@ end subroutine Primitive2ConservativePolyCellsM
!!$ Prim2Con doesn't change for tracers with the addition of a B-field!
!!$
- /*@@
- @routine primitive2conservativegeneralM
- @date Aug 31, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke
- @desc
- Converts primitive to conserved variables for the boundary extended data.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
-
-@@*/
-
-subroutine primitive2conservativegeneralM(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, ierr
- CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl, &
- gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr, &
- vlowx, vlowy, vlowz
-
- ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallPlus)
- ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallMinus)
-
- !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
- !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr)
- do k = GRHydro_stencil , cctk_lsh(3) - GRHydro_stencil + 1
- do j = GRHydro_stencil , cctk_lsh(2) - GRHydro_stencil + 1
- do i = GRHydro_stencil , cctk_lsh(1) - GRHydro_stencil + 1
-
- gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset))
- gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset))
- gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset))
- gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset))
- gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset))
- gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset))
- gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset))
- gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset))
- gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset))
- gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset))
- gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset))
- gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset))
-
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
-
- call prim2conpolytypeM(-1, gxxl,gxyl,gxzl,&
- gyyl,gyzl,gzzl, &
- avg_detl,densminus(i,j,k),sxminus(i,j,k),&
- syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),&
- Bvecxminus(i,j,k),Bvecyminus(i,j,k),Bveczminus(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 prim2conpolytypeM(-1, gxxr,gxyr,gxzr,&
- gyyr,gyzr,gzzr, &
- avg_detr, densplus(i,j,k),sxplus(i,j,k),&
- syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),&
- Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(i,j,k),&
- rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),&
- velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
- w_lorentzplus(i,j,k))
- if (densminus(i,j,k) .ne. densminus(i,j,k)) then
- !$OMP CRITICAL
- call CCTK_WARN(1, "NaN in densminus(i,j,k) (Prim2Con)")
- !$OMP END CRITICAL
- endif
- end do
- end do
- end do
- !$OMP END PARALLEL DO
-
-end subroutine primitive2conservativegeneralM
-
- /*@@
- @routine Primitive2ConservativeCellsGeneralM
- @date Ag 31, 2010
- @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke
- @desc
- Converts primitive to conserved variables for the boundary extended data.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
-
-@@*/
-
-subroutine Primitive2ConservativeCellsGeneralM(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i, j, k, ierr
- CCTK_REAL :: det,gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt, &
- vlowx, vlowy, vlowz
- CCTK_REAL Bdotv,ab0,b2,blowx,blowy,blowz
-
- ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConCellsCall)
-
- 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
-
- gxxpt = gxx(i,j,k)
- gxypt = gxy(i,j,k)
- gxzpt = gxz(i,j,k)
- gyypt = gyy(i,j,k)
- gyzpt = gyz(i,j,k)
- gzzpt = gzz(i,j,k)
-
- det = SPATIAL_DETERMINANT(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt)
-
- w_lorentz(i,j,k)=1.d0/sqrt(1.d0-DOTPT2(velx(i,j,k),vely(i,j,k),velz(i,j,k)))
-
- vlowx = gxxpt * velx(i,j,k) + &
- gxypt * vely(i,j,k) + &
- gxzpt * velz(i,j,k)
- vlowy = gxypt * velx(i,j,k) + &
- gyypt * vely(i,j,k) + &
- gyzpt * velz(i,j,k)
- vlowz = gxzpt * velx(i,j,k) + &
- gyzpt * vely(i,j,k) + &
- gzzpt * velz(i,j,k)
-
- Bdotv=DOTPT(velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k))
- ab0=w_lorentz(i,j,k)*Bdotv
- b2 = DOTPT2(Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k))/w_lorentz(i,j,k)**2+Bdotv**2
-
- blowx = (gxx(i,j,k)*Bvecx(i,j,k) + gxy(i,j,k)*Bvecy(i,j,k) + gxz(i,j,k)*Bvecz(i,j,k))/ &
- w_lorentz(i,j,k) + w_lorentz(i,j,k)*Bdotv*vlowx
- blowy = (gxy(i,j,k)*Bvecx(i,j,k) + gyy(i,j,k)*Bvecy(i,j,k) + gyz(i,j,k)*Bvecz(i,j,k))/ &
- w_lorentz(i,j,k) + w_lorentz(i,j,k)*Bdotv*vlowy
- blowz = (gxz(i,j,k)*Bvecx(i,j,k) + gyz(i,j,k)*Bvecy(i,j,k) + gzz(i,j,k)*Bvecz(i,j,k))/ &
- w_lorentz(i,j,k) + w_lorentz(i,j,k)*Bdotv*vlowz
-
- dens(i,j,k) = sqrt(det) * rho(i,j,k) * w_lorentz(i,j,k)
- sx(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* &
- w_lorentz(i,j,k)*w_lorentz(i,j,k) * vlowx - ab0*blowx)
- sy(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* &
- w_lorentz(i,j,k)*w_lorentz(i,j,k) * vlowy - ab0*blowy)
- sz(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* &
- w_lorentz(i,j,k)*w_lorentz(i,j,k) * vlowz - ab0*blowz)
- tau(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* &
- w_lorentz(i,j,k)*w_lorentz(i,j,k) - press(i,j,k)-b2/2.0-ab0**2) - dens(i,j,k)
-
- end do
- end do
- end do
-
-end subroutine Primitive2ConservativeCellsGeneralM
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90
index 5bce1c7..065d2bd 100644
--- a/src/GRHydro_Reconstruct.F90
+++ b/src/GRHydro_Reconstruct.F90
@@ -875,19 +875,11 @@ subroutine Reconstruction(CCTK_ARGUMENTS)
if (CCTK_EQUALS(recon_vars,"primitive").or.&
CCTK_EQUALS(recon_method,"ppm")) then
- if (use_eosgeneral == 0) then
- if(evolve_mhd.ne.0) then
- call primitive2conservativeM(CCTK_PASS_FTOF)
- else
- call primitive2conservative(CCTK_PASS_FTOF)
- endif
- else
- if(evolve_mhd.ne.0) then
- call primitive2conservativegeneralM(CCTK_PASS_FTOF)
- else
- call primitive2conservativegeneral(CCTK_PASS_FTOF)
- endif
- end if
+ if(evolve_mhd.ne.0) then
+ call primitive2conservativeM(CCTK_PASS_FTOF)
+ else
+ call primitive2conservative(CCTK_PASS_FTOF)
+ endif
else if (CCTK_EQUALS(recon_vars,"conservative")) then
if(evolve_mhd.ne.0) then
call Conservative2PrimitiveBoundsM(CCTK_PASS_FTOF)
diff --git a/src/GRHydro_ReconstructM.F90 b/src/GRHydro_ReconstructM.F90
index 10da759..654d9a7 100644
--- a/src/GRHydro_ReconstructM.F90
+++ b/src/GRHydro_ReconstructM.F90
@@ -675,11 +675,7 @@ subroutine ReconstructionM(CCTK_ARGUMENTS)
endif
if (CCTK_EQUALS(recon_vars,"primitive")) then
- if (use_eosgeneral == 0) then
- call primitive2conservativeM(CCTK_PASS_FTOF)
- else
- call primitive2conservativegeneralM(CCTK_PASS_FTOF)
- end if
+ call primitive2conservativeM(CCTK_PASS_FTOF)
else if (CCTK_EQUALS(recon_vars,"conservative")) then
call Conservative2PrimitiveBoundsM(CCTK_PASS_FTOF)
else
diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90
index b16e483..db78bd9 100644
--- a/src/GRHydro_ReconstructPoly.F90
+++ b/src/GRHydro_ReconstructPoly.F90
@@ -851,19 +851,11 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS)
if (CCTK_EQUALS(recon_vars,"primitive").or.&
CCTK_EQUALS(recon_method,"ppm")) then
- if (use_eosgeneral == 0) then
- if(evolve_mhd.ne.0) then
- call Prim2ConservativePolytypeM(CCTK_PASS_FTOF)
- else
- call Prim2ConservativePolytype(CCTK_PASS_FTOF)
- endif
- else
- if(evolve_mhd.ne.0) then
- call primitive2conservativegeneral(CCTK_PASS_FTOF)
- else
- call primitive2conservativegeneral(CCTK_PASS_FTOF)
- endif
- end if
+ if(evolve_mhd.ne.0) then
+ call Prim2ConservativePolytypeM(CCTK_PASS_FTOF)
+ else
+ call Prim2ConservativePolytype(CCTK_PASS_FTOF)
+ endif
else if (CCTK_EQUALS(recon_vars,"conservative")) then
if(evolve_mhd.ne.0) then
call Con2PrimBoundsPolytype(CCTK_PASS_FTOF)
diff --git a/src/GRHydro_ReconstructPolyM.F90 b/src/GRHydro_ReconstructPolyM.F90
index 88ef183..70d2b12 100644
--- a/src/GRHydro_ReconstructPolyM.F90
+++ b/src/GRHydro_ReconstructPolyM.F90
@@ -652,14 +652,9 @@ subroutine ReconstructionPolytypeM(CCTK_ARGUMENTS)
endif
if (CCTK_EQUALS(recon_vars,"primitive")) then
-
- if (use_eosgeneral == 0) then
- call Prim2ConservativePolytypeM(CCTK_PASS_FTOF)
- else
- call primitive2conservativegeneralM(CCTK_PASS_FTOF)
- end if
+ call Prim2ConservativePolytypeM(CCTK_PASS_FTOF)
else if (CCTK_EQUALS(recon_vars,"conservative")) then
- call Con2PrimBoundsPolytypeM(CCTK_PASS_FTOF)
+ call Con2PrimBoundsPolytypeM(CCTK_PASS_FTOF)
else
call CCTK_WARN(0,"Variable type to reconstruct not recognized.")
end if
diff --git a/src/GRHydro_RegisterGZ.cc b/src/GRHydro_RegisterGZ.cc
deleted file mode 100644
index 6029145..0000000
--- a/src/GRHydro_RegisterGZ.cc
+++ /dev/null
@@ -1,182 +0,0 @@
-// register.cc -- register variables with various thorns that need-to-know
-// $Header$
-//
-// GRHydro_register_GZPatchSystem - register with GZPatchSystem
-//
-// Cut 'n paste job from BackgroundWaveToy, author J Thornburg...
-
-#include <cstdio>
-#include <string>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-using namespace std;
-
-//******************************************************************************
-
-//
-// This function is called by the Cactus scheduler to (maybe) register
-// our to-be-interpatch-synchronized variables with GZPatchSystem. It
-// checks if the GZPatchSystem registration function has been provided
-// (which it should be if and only if GZPatchSystem is active), and if so,
-// does the registration.
-//
-// If we're using Carpet, this function must be called in meta mode.
-//
-extern "C"void GRHydro_register_GZPatchSystem(CCTK_ARGUMENTS)
-{
- DECLARE_CCTK_PARAMETERS;
-
- if (CCTK_IsFunctionAliased("GZPatchSystem_register_sync"))
- {
- CCTK_VInfo(CCTK_THORNSTRING,
- "registering to-be-interpatch-synchronized variables "
- "with GZPatchSystem");
-
- if(CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) {
- if(clean_divergence) {
- string var[10] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec", "GRHydro::psidc"};
- for (int i = 0; i < 10; i++)
- {
- int status =
- GZPatchSystem_register_sync(var[i].c_str());
- if (status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-synchronized!\n"
- " (GZPatchSystem_register_sync() error code %d)\n",
- var[i].c_str(), status);
- }
- }
- } else {
- string var[9] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec"};
- for (int i = 0; i < 9; i++)
- {
- int status =
- GZPatchSystem_register_sync(var[i].c_str());
- if (status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-synchronized!\n"
- " (GZPatchSystem_register_sync() error code %d)\n",
- var[i].c_str(), status);
- }
- }
- }
- } else {
- string var[8] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon"};
- for (int i = 0; i < 8; i++)
- {
- int status =
- GZPatchSystem_register_sync(var[i].c_str());
- if (status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-synchronized!\n"
- " (GZPatchSystem_register_sync() error code %d)\n",
- var[i].c_str(), status);
- }
- }
- }
- }
- else
- {
- CCTK_WARN(1, "Function GZPatchSystem_register_sync not registered!");
- }
-
- if (CCTK_IsFunctionAliased("GZPatchSystem_register_cxform"))
- {
- CCTK_VInfo(CCTK_THORNSTRING,
- "registering to-be-cxformed variables with GZPatchSystem");
-
- if(CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) {
- if(clean_divergence) {
- string var[13] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec", "GRHydro::psidc",
- "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"};
- for (int i = 0; i < 13; i++)
- for (int j = 0; j < 3; j++)
- {
- int ps_status =
- GZPatchSystem_register_cxform(j, var[i].c_str());
- if (ps_status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-cxformhronized!\n"
- " (GZPatchSystem_register_cxform() error code %d)\n",
- var[i].c_str(), ps_status);
- }
- }
- } else {
- string var[12] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec",
- "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"};
- for (int i = 0; i < 12; i++)
- for (int j = 0; j < 3; j++)
- {
- int ps_status =
- GZPatchSystem_register_cxform(j, var[i].c_str());
- if (ps_status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-cxformhronized!\n"
- " (GZPatchSystem_register_cxform() error code %d)\n",
- var[i].c_str(), ps_status);
- }
- }
- }
-
- } else {
- string var[11] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon",
- "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"};
- for (int i = 0; i < 11; i++)
- for (int j = 0; j < 3; j++)
- {
- int ps_status =
- GZPatchSystem_register_cxform(j, var[i].c_str());
- if (ps_status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-cxformhronized!\n"
- " (GZPatchSystem_register_cxform() error code %d)\n",
- var[i].c_str(), ps_status);
- }
- }
- }
- }
- else
- {
- CCTK_WARN(1, "Function GZPatchSystem_register_cxform not registered!");
- }
-
-}
diff --git a/src/GRHydro_RegisterGZM.cc b/src/GRHydro_RegisterGZM.cc
deleted file mode 100644
index 2e8ab81..0000000
--- a/src/GRHydro_RegisterGZM.cc
+++ /dev/null
@@ -1,138 +0,0 @@
-// register.cc -- register variables with various thorns that need-to-know
-// $Header$
-//
-// GRHydro_register_GZPatchSystem - register with GZPatchSystem
-//
-// Cut 'n paste job from BackgroundWaveToy, author J Thornburg...
-
-#include <cstdio>
-#include <string>
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-using namespace std;
-
-//******************************************************************************
-
-//
-// This function is called by the Cactus scheduler to (maybe) register
-// our to-be-interpatch-synchronized variables with GZPatchSystem. It
-// checks if the GZPatchSystem registration function has been provided
-// (which it should be if and only if GZPatchSystem is active), and if so,
-// does the registration.
-//
-// If we're using Carpet, this function must be called in meta mode.
-//
-extern "C"void GRHydro_register_GZPatchSystemM(CCTK_ARGUMENTS)
-{
- DECLARE_CCTK_PARAMETERS;
-
- if (CCTK_IsFunctionAliased("GZPatchSystem_register_sync"))
- {
- CCTK_VInfo(CCTK_THORNSTRING,
- "registering to-be-interpatch-synchronized variables "
- "with GZPatchSystem");
-
- if(clean_divergence) {
- string var[10] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec", "GRHydro::psidc"};
- for (int i = 0; i < 10; i++)
- {
- int status =
- GZPatchSystem_register_sync(var[i].c_str());
- if (status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-synchronized!\n"
- " (GZPatchSystem_register_sync() error code %d)\n",
- var[i].c_str(), status);
- }
- }
- } else {
- string var[9] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec"};
- for (int i = 0; i < 9; i++)
- {
- int status =
- GZPatchSystem_register_sync(var[i].c_str());
- if (status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-synchronized!\n"
- " (GZPatchSystem_register_sync() error code %d)\n",
- var[i].c_str(), status);
- }
- }
-
- }
- }
- else
- {
- CCTK_WARN(1, "Function GZPatchSystem_register_sync not registered!");
- }
-
- if (CCTK_IsFunctionAliased("GZPatchSystem_register_cxform"))
- {
- CCTK_VInfo(CCTK_THORNSTRING,
- "registering to-be-cxformed variables with GZPatchSystem");
-
- if(clean_divergence) {
- string var[13] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec", "GRHydro::psidc",
- "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"};
- for (int i = 0; i < 13; i++)
- for (int j = 0; j < 3; j++)
- {
- int ps_status =
- GZPatchSystem_register_cxform(j, var[i].c_str());
- if (ps_status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-cxformhronized!\n"
- " (GZPatchSystem_register_cxform() error code %d)\n",
- var[i].c_str(), ps_status);
- }
- }
- } else {
- string var[12] = {"HydroBase::rho", "HydroBase::press", "HydroBase::eps",
- "HydroBase::vel",
- "GRHydro::dens", "GRHydro::tau", "HydroBase::w_lorentz",
- "GRHydro::scon", "HydroBase::Bvec",
- "ADMBase::metric", "ADMBase::curv", "ADMBase::shift"};
- for (int i = 0; i < 12; i++)
- for (int j = 0; j < 3; j++)
- {
- int ps_status =
- GZPatchSystem_register_cxform(j, var[i].c_str());
- if (ps_status < 0)
- {
- CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "***** GRHydro_register_GZPatchSystem():\n"
- " error registering var group %s to be "
- "interpatch-cxformhronized!\n"
- " (GZPatchSystem_register_cxform() error code %d)\n",
- var[i].c_str(), ps_status);
- }
- }
- }
- }
- else
- {
- CCTK_WARN(1, "Function GZPatchSystem_register_cxform not registered!");
- }
-
-}
diff --git a/src/GRHydro_RiemannSolve.F90 b/src/GRHydro_RiemannSolve.F90
index 84c01e2..be2c2d4 100644
--- a/src/GRHydro_RiemannSolve.F90
+++ b/src/GRHydro_RiemannSolve.F90
@@ -127,235 +127,4 @@ end subroutine RiemannSolvePolytype
- /*@@
- @routine RiemannSolveGeneral
- @date Tue Mar 19 11:40:20 2002
- @author Ian Hawke
- @desc
- The Riemann solvers for the new general EOS routines.
- This sets the fluxes from the left and right reconstructed
- states, so that after this routine they are effectively
- scratch space.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
-
-@@*/
-
-
-subroutine RiemannSolveGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i,j,k, ierr
- CCTK_REAL, dimension(5) :: tmp_flux, cons_p, cons_m
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL :: avg_det, avg_alp, avg_beta
- CCTK_REAL :: psi4h, gxxh, gyyh, gzzh, gxyh, gxzh, gyzh
-
- densflux = 0.d0
- sxflux = 0.d0
- syflux = 0.d0
- szflux = 0.d0
- tauflux = 0.d0
-
-!!$ Do the EOS call to set the pressure, derivative and cs2
-
- ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallPlus)
- ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallMinus)
-
-
- if (CCTK_EQUALS(riemann_solver,"HLLE")) then
-
- call GRHydro_HLLEGeneral(CCTK_PASS_FTOF)
-
- if (evolve_tracer .ne. 0) then
-
- call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF)
-
- end if
-
- else
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
-!!$ Set the left (p for plus) and right (m for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
- prim_p(6) = pressplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
-
-!!$ Set metric terms at interface
-
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + &
- gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + &
- gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + &
- gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + &
- gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + &
- gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + &
- gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
- if (flux_direction == 1) then
-
- call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
- tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), &
- prim_p(2),prim_p(6), &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = 0.5d0 * tmp_flux(5)
-
- call num_x_flux(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
- tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), &
- prim_m(2),prim_m(6), &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5)
-
- else if (flux_direction == 2) then
-
- call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
- tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5), &
- prim_p(3),prim_p(6), &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = 0.5d0 * tmp_flux(5)
-
- call num_x_flux(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
- tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5), &
- prim_m(3),prim_m(6), &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5)
-
- else if (flux_direction == 3) then
-
- call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
- tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
- prim_p(4),prim_p(6), &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = 0.5d0 * tmp_flux(5)
-
- call num_x_flux(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
- tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
- prim_m(4),prim_m(6), &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5)
-
- else
-
- call CCTK_WARN(0, "Flux direction not x,y,z")
-
- end if
-
- end do
- end do
- end do
-
- if (CCTK_EQUALS(riemann_solver,"Roe")) then
-
- call GRHydro_RoeSolveGeneral(CCTK_PASS_FTOF)
-
- if (evolve_tracer .ne. 0) then
-
- call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF)
-
- end if
-
- else if (CCTK_EQUALS(riemann_solver,"Marquina")) then
-
- call GRHydro_MarquinaGeneral(CCTK_PASS_FTOF)
-
-!!$ Tracers are built directly in to the Marquina solver
-
- end if
-
- end if
-
-end subroutine RiemannSolveGeneral
-
diff --git a/src/GRHydro_RiemannSolveM.F90 b/src/GRHydro_RiemannSolveM.F90
index 8944b26..bd16778 100644
--- a/src/GRHydro_RiemannSolveM.F90
+++ b/src/GRHydro_RiemannSolveM.F90
@@ -140,358 +140,5 @@ end subroutine RiemannSolvePolytypeM
-/*@@
-@routine RiemannSolveGeneralM
-@date Tue Mar 19 11:40:20 2002
-@author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
-@desc
-The Riemann solvers for the new general EOS routines.
-This sets the fluxes from the left and right reconstructed
-states, so that after this routine they are effectively
-scratch space.
-@enddesc
-@calls
-@calledby
-@history
-
-@endhistory
-
-@@*/
-
-
-subroutine RiemannSolveGeneralM(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_INT :: i,j,k, ierr
- CCTK_REAL, dimension(8) :: tmp_flux, cons_p, cons_m
- CCTK_REAL, dimension(6) :: prim_p, prim_m
- CCTK_REAL :: avg_det, avg_alp, avg_beta
- CCTK_REAL :: gxxh, gyyh, gzzh, gxyh, gxzh, gyzh
- CCTK_REAL :: avg_betax, avg_betay, avg_betaz
- CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm
- CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p
- CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm
- CCTK_REAL :: pressstarp,pressstarm,rhoenth_p,rhoenth_m
- CCTK_REAL :: velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm
-
- CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm
-
- densflux = 0.d0
- sxflux = 0.d0
- syflux = 0.d0
- szflux = 0.d0
- tauflux = 0.d0
- Bvecxflux = 0.d0
- Bvecyflux = 0.d0
- Bveczflux = 0.d0
- if(clean_divergence.ne.0) then
- psidcflux = 0.d0
- endif
-
-!!$ Do the EOS call to set the pressure, derivative and cs2
-
- ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallPlus)
- ierr = EOS_SetGFs(cctkGH, EOS_RiemannCallMinus)
-
-
- if (CCTK_EQUALS(riemann_solver,"HLLE")) then
-
- call GRHydro_HLLEGeneralM(CCTK_PASS_FTOF)
-
- if (evolve_tracer .ne. 0) then
-
-!!$ No b-field component for tracers!
- call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF)
-
- end if
-
- else
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
-!!$ Set the left (p for plus) and right (m for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
- cons_p(6) = Bvecxplus(i,j,k)
- cons_p(7) = Bvecyplus(i,j,k)
- cons_p(8) = Bveczplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(6) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(7) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(8) = Bveczminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
- prim_p(6) = pressplus(i,j,k)
- rhoenth_p = prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
- rhoenth_m = prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)
-
- if(clean_divergence.ne.0) then
- psidcp = psidcplus(i,j,k)
- psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset)
- endif
-
-!!$ Set metric terms at interface
-
- if (shift_state .ne. 0) then
- avg_betax = 0.5d0*(betax(i+xoffset,j+yoffset,k+zoffset)+betax(i,j,k))
- avg_betay = 0.5d0*(betay(i+xoffset,j+yoffset,k+zoffset)+betay(i,j,k))
- avg_betaz = 0.5d0*(betaz(i+xoffset,j+yoffset,k+zoffset)+betaz(i,j,k))
- if (flux_direction == 1) then
- avg_beta = avg_betax
- else if (flux_direction == 2) then
- avg_beta = avg_betay
- else if (flux_direction == 3) then
- avg_beta = avg_betaz
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- avg_betax = 0.d0
- avg_betay = 0.d0
- avg_betaz = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + &
- gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + &
- gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + &
- gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + &
- gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + &
- gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + &
- gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
- vxtp = prim_p(2)-avg_betax/avg_alp
- vytp = prim_p(3)-avg_betay/avg_alp
- vztp = prim_p(4)-avg_betaz/avg_alp
- vxtm = prim_m(2)-avg_betax/avg_alp
- vytm = prim_m(3)-avg_betay/avg_alp
- vztm = prim_m(4)-avg_betaz/avg_alp
-
- call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- prim_p(2),prim_p(3),prim_p(4),cons_p(6),cons_p(7),cons_p(8), &
- velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, &
- bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp)
- call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- prim_m(2),prim_m(3),prim_m(4),cons_m(6),cons_m(7),cons_m(8), &
- velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, &
- bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm)
-
- ab0p = wp*bdotvp
- ab0m = wm*bdotvm
-
- vA2p = b2p/(rhoenth_p+b2p)
- vA2m = b2m/(rhoenth_m+b2m)
-
-!!$ p^* = p+b^2/2 in Anton et al.
- pressstarp = prim_p(6)+0.5d0*b2p
- pressstarm = prim_m(6)+0.5d0*b2m
-
- if (flux_direction == 1) then
-
- call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),&
- cons_p(6),cons_p(7),cons_p(8),&
- tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5), &
- tmp_flux(6),tmp_flux(7),tmp_flux(8), &
- vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = 0.5d0 * tmp_flux(5)
- Bvecxflux(i,j,k) = 0.5d0 * tmp_flux(6)
- Bvecyflux(i,j,k) = 0.5d0 * tmp_flux(7)
- Bveczflux(i,j,k) = 0.5d0 * tmp_flux(8)
-
- if(clean_divergence.ne.0) then
- psidcf = cons_p(6)
- psidcflux(i,j,k) = 0.5d0 * psidcf
- endif
-
- call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),&
- cons_m(6),cons_m(7),cons_m(8),&
- tmp_flux(1),tmp_flux(2),tmp_flux(3),tmp_flux(4),tmp_flux(5),&
- tmp_flux(6),tmp_flux(7),tmp_flux(8),&
- vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5)
- Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
- Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
- Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-
- if(clean_divergence.ne.0) then
- psidcf = cons_m(6)
- psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
- endif
-
- else if (flux_direction == 2) then
-
- call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),&
- cons_p(7),cons_p(8),cons_p(6),&
- tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),&
- tmp_flux(7),tmp_flux(8),tmp_flux(6),&
- vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = 0.5d0 * tmp_flux(5)
- Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
- Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
- Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-
- if(clean_divergence.ne.0) then
- psidcf = cons_p(7)
- psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
- endif
-
- call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),&
- cons_m(7),cons_m(8),cons_m(6),&
- tmp_flux(1),tmp_flux(3),tmp_flux(4),tmp_flux(2),tmp_flux(5),&
- tmp_flux(7),tmp_flux(8),tmp_flux(6),&
- vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5)
- Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
- Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
- Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-
- if(clean_divergence.ne.0) then
- psidcf = cons_m(7)
- psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
- endif
-
- else if (flux_direction == 3) then
-
- call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),&
- cons_p(8),cons_p(6),cons_p(7),&
- tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5),&
- tmp_flux(8),tmp_flux(6),tmp_flux(7), &
- vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = 0.5d0 * tmp_flux(5)
- Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
- Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
- Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-
- if(clean_divergence.ne.0) then
- psidcf = cons_p(8)
- psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
- endif
-
- call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),&
- cons_m(8),cons_m(6),cons_m(7),&
- tmp_flux(1),tmp_flux(4),tmp_flux(2),tmp_flux(3),tmp_flux(5), &
- tmp_flux(8),tmp_flux(6),tmp_flux(7), &
- vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
- avg_det,avg_alp,avg_beta)
-
- densflux(i,j,k) = densflux(i,j,k) + 0.5d0 * tmp_flux(1)
- sxflux(i,j,k) = sxflux(i,j,k) + 0.5d0 * tmp_flux(2)
- syflux(i,j,k) = syflux(i,j,k) + 0.5d0 * tmp_flux(3)
- szflux(i,j,k) = szflux(i,j,k) + 0.5d0 * tmp_flux(4)
- tauflux(i,j,k) = tauflux(i,j,k) + 0.5d0 * tmp_flux(5)
- Bvecxflux(i,j,k)= Bvecxflux(i,j,k)+ 0.5d0 * tmp_flux(6)
- Bvecyflux(i,j,k)= Bvecyflux(i,j,k)+ 0.5d0 * tmp_flux(7)
- Bveczflux(i,j,k)= Bveczflux(i,j,k)+ 0.5d0 * tmp_flux(8)
-
- if(clean_divergence.ne.0) then
- psidcf = cons_m(8)
- psidcflux(i,j,k) = psidcflux(i,j,k) + 0.5d0 * psidcf
- endif
-
- else
-
- call CCTK_WARN(0, "Flux direction not x,y,z")
-
- end if
-
- end do
- end do
- end do
-
- if (CCTK_EQUALS(riemann_solver,"Roe")) then
-
- call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!")
-
-!!$
-!!$ call GRHydro_RoeSolveGeneral(CCTK_PASS_FTOF)
-!!$
-!!$ if (evolve_tracer .ne. 0) then
-!!$
-!!$ call GRHydro_HLLE_TracerGeneral(CCTK_PASS_FTOF)
-!!$
-!!$ end if
-!!$
- else if (CCTK_EQUALS(riemann_solver,"Marquina")) then
-
- call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!")
-
-!!$
-!!$ call GRHydro_MarquinaGeneral(CCTK_PASS_FTOF)
-!!$
-!!$ Tracers are built directly in to the Marquina solver
-
- end if
-
- end if
-
-end subroutine RiemannSolveGeneralM
diff --git a/src/GRHydro_RoeSolver.F90 b/src/GRHydro_RoeSolver.F90
index 329c081..26f9999 100644
--- a/src/GRHydro_RoeSolver.F90
+++ b/src/GRHydro_RoeSolver.F90
@@ -310,246 +310,3 @@ subroutine GRHydro_RoeSolve(CCTK_ARGUMENTS)
end subroutine GRHydro_RoeSolve
- /*@@
- @routine GRHydro_RoeSolveGeneral
- @date Sat Jan 26 01:55:55 2002
- @author Pedro Montero, Ian Hawke
- @desc
- Wrapper routine to calculate the Roe fluxes and hence the update
- terms.
- @enddesc
- @calls
- @calledby
- @history
-
- @endhistory
-
-@@*/
-
-subroutine GRHydro_RoeSolveGeneral(CCTK_ARGUMENTS)
-
- USE GRHydro_Scalars
- use GRHydro_Eigenproblem
-
- implicit none
-
- DECLARE_CCTK_ARGUMENTS
- DECLARE_CCTK_PARAMETERS
- DECLARE_CCTK_FUNCTIONS
-
- CCTK_REAL, dimension(5) :: roeflux,roeave,qdiff,cons_p,cons_m,&
- cons_ave
- CCTK_REAL, dimension(6) :: prim_p, prim_m, prim_ave
- CCTK_REAL :: avg_alp,avg_beta,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- avg_det,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,&
- usendh,alp_l,alp_r,psi4h,cs2_ave,dpdeps_ave
- CCTK_INT :: i,j,k,m
-
- CCTK_INT :: type_bits, trivial, not_trivial, ierr
-
- if (flux_direction == 1) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", &
- &"trivial")
- else if (flux_direction == 2) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", &
- &"trivial")
- else if (flux_direction == 3) then
- call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ")
- call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", &
- &"trivial")
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(1) = rhoplus(i,j,k)
- prim_p(2) = velxplus(i,j,k)
- prim_p(3) = velyplus(i,j,k)
- prim_p(4) = velzplus(i,j,k)
- prim_p(5) = epsplus(i,j,k)
-
- prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset)
- prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_p(6) = pressplus(i,j,k)
-
- prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset)
-
-!!$ Set the Roe average of the fluid variables
-
- call roeaverage(prim_p(1:5), prim_m(1:5), roeave)
-
- rho_ave(i,j,k) = roeave(1)
- velx_ave(i,j,k) = roeave(2)
- vely_ave(i,j,k) = roeave(3)
- velz_ave(i,j,k) = roeave(4)
- eps_ave(i,j,k) = roeave(5)
-
- end do
- end do
- end do
-
- ierr = EOS_SetGFs(cctkGH, EOS_RoeAverageCall)
-
- do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil
-
-!!$ Set the left (p for plus) and right (m_i for minus, i+1) states
-
- cons_p(1) = densplus(i,j,k)
- cons_p(2) = sxplus(i,j,k)
- cons_p(3) = syplus(i,j,k)
- cons_p(4) = szplus(i,j,k)
- cons_p(5) = tauplus(i,j,k)
-
- cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset)
- cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset)
-
- prim_ave(1) = rho_ave(i,j,k)
- prim_ave(2) = velx_ave(i,j,k)
- prim_ave(3) = vely_ave(i,j,k)
- prim_ave(4) = velz_ave(i,j,k)
- prim_ave(5) = eps_ave(i,j,k)
- prim_ave(6) = press_ave(i,j,k)
- cs2_ave = eos_cs2_ave(i,j,k)
- dpdeps_ave = eos_dpdeps_ave(i,j,k)
-
- roeflux = 0.d0
- qdiff = 0.d0
-
-!!$ Calculate jumps in conserved variables
-
- do m = 1,5
- qdiff(m) = cons_m(m) - cons_p(m)
- end do
-
-!!$ Set metric terms at interface
-
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
- else
- avg_beta = 0.d0
- end if
-
- avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
-
- gxxh = 0.5d0 * (gxx(i+xoffset,j+yoffset,k+zoffset) + &
- gxx(i,j,k))
- gxyh = 0.5d0 * (gxy(i+xoffset,j+yoffset,k+zoffset) + &
- gxy(i,j,k))
- gxzh = 0.5d0 * (gxz(i+xoffset,j+yoffset,k+zoffset) + &
- gxz(i,j,k))
- gyyh = 0.5d0 * (gyy(i+xoffset,j+yoffset,k+zoffset) + &
- gyy(i,j,k))
- gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + &
- gyz(i,j,k))
- gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + &
- gzz(i,j,k))
-
- avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
-
-!!$ If the Riemann problem is trivial the flux is already correct
-
- if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then
-
- else !!! The end of this branch is right at the bottom of the routine
-
- call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
- avg_det,gxxh, gxyh, gxzh, gyyh, gyzh, gzzh)
-
- if (flux_direction == 1) then
- usendh = uxxh
- else if (flux_direction == 2) then
- usendh = uyyh
- else if (flux_direction == 3) then
- usendh = uzzh
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
-!!$ Convert to conserved variables and find the part of the Roe
-!!$ flux that requires the spectral decomposition.
-!!$ The conversion to conserved variables is just to set the
-!!$ pressure at this point (means this routine doesn''t need
-!!$ the EOS interface).
-
- if (flux_direction == 1) then
- call eigenproblem_general(&
- prim_ave(1),prim_ave(2),prim_ave(3),prim_ave(4),prim_ave(5), &
- prim_ave(6),cs2_ave,dpdeps_ave, &
- gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
- usendh,avg_alp,avg_beta, &
- qdiff(1),qdiff(2),qdiff(3),qdiff(4),qdiff(5), &
- roeflux(1),roeflux(2),roeflux(3),roeflux(4),roeflux(5))
- else if (flux_direction == 2) then
- call eigenproblem_general(&
- prim_ave(1),prim_ave(3),prim_ave(4),prim_ave(2),prim_ave(5), &
- prim_ave(6),cs2_ave,dpdeps_ave, &
- gyyh,gyzh,gxyh,gzzh,gxzh,gxxh, &
- usendh,avg_alp,avg_beta, &
- qdiff(1),qdiff(3),qdiff(4),qdiff(2),qdiff(5), &
- roeflux(1),roeflux(3),roeflux(4),roeflux(2),roeflux(5))
- else if (flux_direction == 3) then
- call eigenproblem_general(&
- prim_ave(1),prim_ave(4),prim_ave(2),prim_ave(3),prim_ave(5), &
- prim_ave(6),cs2_ave,dpdeps_ave, &
- gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, &
- usendh,avg_alp,avg_beta, &
- qdiff(1),qdiff(4),qdiff(2),qdiff(3),qdiff(5), &
- roeflux(1),roeflux(4),roeflux(2),roeflux(3),roeflux(5))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
-
-!!$ Roe flux
-
- densflux(i,j,k) = densflux(i,j,k) - 0.5d0 * roeflux(1)
- sxflux(i,j,k) = sxflux(i,j,k) - 0.5d0 * roeflux(2)
- syflux(i,j,k) = syflux(i,j,k) - 0.5d0 * roeflux(3)
- szflux(i,j,k) = szflux(i,j,k) - 0.5d0 * roeflux(4)
- tauflux(i,j,k) = tauflux(i,j,k) - 0.5d0 * roeflux(5)
-
- end if !!! The end of the SpaceMask check for a trivial RP.
-
- enddo
- enddo
- enddo
-
-end subroutine GRHydro_RoeSolveGeneral
diff --git a/src/GRHydro_Scalars.F90 b/src/GRHydro_Scalars.F90
index 2cc245e..ae03e02 100644
--- a/src/GRHydro_Scalars.F90
+++ b/src/GRHydro_Scalars.F90
@@ -20,11 +20,5 @@
LOGICAL :: ANALYTICAL
LOGICAL :: FAST
- CCTK_INT :: EOS_RiemannCallPlus, EOS_RiemannCallMinus, EOS_RoeAverageCall
- CCTK_INT :: EOS_Con2PrimCall
- CCTK_INT :: EOS_Prim2ConBndCallPlus, EOS_Prim2ConBndCallMinus
- CCTK_INT :: EOS_Prim2ConCellsCall
-
- CCTK_REAL :: eosgeneral_pmin
end module GRHydro_Scalars
diff --git a/src/GRHydro_Startup.F90 b/src/GRHydro_Startup.F90
index 444e944..5b6c97c 100644
--- a/src/GRHydro_Startup.F90
+++ b/src/GRHydro_Startup.F90
@@ -38,391 +38,10 @@ integer function GRHydro_Startup()
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
- CCTK_INT :: ierr, flags, table_handle, var, strlen, strend, n_base_indep_gfs
- CCTK_INT, dimension(:), allocatable :: indep_gfs
- CCTK_INT, dimension(4) :: dep_gfs
- character(len=256) :: eos_name_param, var_names_param, var_names, var_gfs
+ CCTK_INT :: ierr
call CCTK_RegisterBanner(ierr, "GRHydro: relativistic hydrodynamics, no ice.")
-!!$ The number of independent variables is different in the polytropic or general case
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- n_base_indep_gfs = 1
- else
- n_base_indep_gfs = 2
- end if
-
- allocate(indep_gfs(n_base_indep_gfs + eosgeneral_n_indeps))
-
-
-!!$ Set up EOS calls.
-
- if (use_eosgeneral .ne. 0) then
-
-!!$ First, set up the minimum pressure
-
- eosgeneral_pmin = 1.d-28
-
-!!$ First, "Plus" boundary extended
-
- flags = UTIL_TABLE_FLAGS_CASE_INSENSITIVE
-
- call Util_TableCreate(table_handle, flags)
-
- call CCTK_FortranString(strlen, eosgeneral_name, eos_name_param)
- eos_name_param = trim(eos_name_param)
- call Util_TableSetString(ierr, table_handle, &
- eos_name_param, "EOS Name")
- call CCTK_FortranString(strlen, eosgeneral_indep_names, &
- var_names_param)
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- var_names = "Rho "//var_names_param(1:strlen)
- else
- var_names = "Rho SpecificInternalEnergy "//var_names_param(1:strlen)
- end if
- var_names = trim(var_names)
- call Util_TableSetString(ierr, table_handle, &
- var_names, "Independent variable names")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure DPressureDSpecificInternalEnergy c_s^2", &
- "Dependent variable names")
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetInt(ierr, table_handle, &
- 1 + eosgeneral_n_indeps, "N independent variables")
- else
- call Util_TableSetInt(ierr, table_handle, &
- 2 + eosgeneral_n_indeps, "N independent variables")
- end if
- call Util_TableSetInt(ierr, table_handle, &
- 3, "N dependent variables")
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhoplus")
- else
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhoplus")
- call CCTK_VarIndex(indep_gfs(2), "GRHydro::epsplus")
- end if
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::pressplus")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::eos_dpdeps_p")
- call CCTK_VarIndex(dep_gfs(3), "GRHydro::eos_cs2_p")
-
- call CCTK_FortranString(strlen, eosgeneral_indep_plus_gfs, var_gfs)
- var_gfs = trim(var_gfs)
-
- strend = 1
-
- do var = 1, eosgeneral_n_indeps
-
- strend = scan(var_gfs, " ")
- if (strend > 0) then
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs(1:strend) )
- var_gfs = var_gfs(strend+1:)
- else
- if (var .ne. eosgeneral_n_indeps) then
- call CCTK_WARN(0, "Something wrong in the eos string")
- end if
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs )
- end if
-
- end do
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- end if
- call Util_TableSetIntArray(ierr, table_handle, 3, dep_gfs, &
- "Dependent GFs")
-
- EOS_RiemannCallPlus = EOS_SetupCall(table_handle)
-
-!!$ Second, "Minus" boundary extended
-!!$ Many of the table entries are the same
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhominus")
- else
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhominus")
- call CCTK_VarIndex(indep_gfs(2), "GRHydro::epsminus")
- end if
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::pressminus")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::eos_dpdeps_m")
- call CCTK_VarIndex(dep_gfs(3), "GRHydro::eos_cs2_m")
-
- call CCTK_FortranString(strlen, eosgeneral_indep_minus_gfs, var_gfs)
- var_gfs = trim(var_gfs)
-
- strend = 1
-
- do var = 1, eosgeneral_n_indeps
-
- strend = scan(var_gfs, " ")
- if (strend > 0) then
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs(1:strend) )
- var_gfs = var_gfs(strend+1:)
- else
- if (var .ne. eosgeneral_n_indeps) then
- call CCTK_WARN(0, "Something wrong in the eos string")
- end if
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs )
- end if
-
- end do
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- end if
- call Util_TableSetIntArray(ierr, table_handle, 3, dep_gfs, &
- "Dependent GFs")
-
- EOS_RiemannCallMinus = EOS_SetupCall(table_handle)
-
-!!$ Third, the "Average" state for the Roe solver
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rho_ave")
- else
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rho_ave")
- call CCTK_VarIndex(indep_gfs(2), "GRHydro::eps_ave")
- end if
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::press_ave")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::eos_dpdeps_ave")
- call CCTK_VarIndex(dep_gfs(3), "GRHydro::eos_cs2_ave")
-
- call CCTK_FortranString(strlen, eosgeneral_indep_gfs, var_gfs)
- var_gfs = trim(var_gfs)
-
- strend = 1
-
- do var = 1, eosgeneral_n_indeps
-
- strend = scan(var_gfs, " ")
- if (strend > 0) then
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs(1:strend) )
- var_gfs = var_gfs(strend+1:)
- else
- if (var .ne. eosgeneral_n_indeps) then
- call CCTK_WARN(0, "Something wrong in the eos string")
- end if
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs )
- end if
-
- end do
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- end if
- call Util_TableSetIntArray(ierr, table_handle, 3, dep_gfs, &
- "Dependent GFs")
-
- EOS_RoeAverageCall = EOS_SetupCall(table_handle)
-
-!!$ Fourth, the "Con2Prim" call.
-!!$ This sets different variables.
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "HydroBase::rho")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure DPressureDSpecificInternalEnergy DPressureDRho SpecificInternalEnergy", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "HydroBase::press")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::eos_dpdeps_temp")
- call CCTK_VarIndex(dep_gfs(3), "GRHydro::eos_dpdrho_temp")
- call CCTK_VarIndex(dep_gfs(4), "HydroBase::eps")
- else
- call CCTK_VarIndex(indep_gfs(1), "HydroBase::rho")
- call CCTK_VarIndex(indep_gfs(2), "HydroBase::eps")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure DPressureDSpecificInternalEnergy DPressureDRho", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "HydroBase::press")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::eos_dpdeps_temp")
- call CCTK_VarIndex(dep_gfs(3), "GRHydro::eos_dpdrho_temp")
- end if
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 4, dep_gfs, &
- "Dependent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 3, dep_gfs, &
- "Dependent GFs")
- end if
-
- EOS_Con2PrimCall = EOS_SetupCall(table_handle)
-
-!!$ Fifth, the "Prim2Con" (boundary) calls.
-!!$ Again, different variables are set.
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhominus")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure SpecificInternalEnergy", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::pressminus")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::epsminus")
- else
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhominus")
- call CCTK_VarIndex(indep_gfs(2), "GRHydro::epsminus")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::pressminus")
- end if
-
- call CCTK_FortranString(strlen, eosgeneral_indep_minus_gfs, var_gfs)
- var_gfs = trim(var_gfs)
-
- strend = 1
-
- do var = 1, eosgeneral_n_indeps
-
- strend = scan(var_gfs, " ")
- if (strend > 0) then
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs(1:strend) )
- var_gfs = var_gfs(strend+1:)
- else
- if (var .ne. eosgeneral_n_indeps) then
- call CCTK_WARN(0, "Something wrong in the eos string")
- end if
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs )
- end if
-
- end do
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 2, dep_gfs, &
- "Dependent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 1, dep_gfs, &
- "Dependent GFs")
- end if
-
- EOS_Prim2ConBndCallMinus = EOS_SetupCall(table_handle)
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhoplus")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure SpecificInternalEnergy", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::pressplus")
- call CCTK_VarIndex(dep_gfs(2), "GRHydro::epsplus")
- else
- call CCTK_VarIndex(indep_gfs(1), "GRHydro::rhoplus")
- call CCTK_VarIndex(indep_gfs(2), "GRHydro::epsplus")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "GRHydro::pressplus")
- end if
-
- call CCTK_FortranString(strlen, eosgeneral_indep_plus_gfs, var_gfs)
- var_gfs = trim(var_gfs)
-
- strend = 1
-
- do var = 1, eosgeneral_n_indeps
-
- strend = scan(var_gfs, " ")
- if (strend > 0) then
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs(1:strend) )
- var_gfs = var_gfs(strend+1:)
- else
- if (var .ne. eosgeneral_n_indeps) then
- call CCTK_WARN(0, "Something wrong in the eos string")
- end if
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs )
- end if
-
- end do
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 2, dep_gfs, &
- "Dependent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 1, dep_gfs, &
- "Dependent GFs")
- end if
-
- EOS_Prim2ConBndCallPlus = EOS_SetupCall(table_handle)
-
-!!$ Sixth, the Prim2Con (cell centre) calls
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call CCTK_VarIndex(indep_gfs(1), "HydroBase::rho")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure SpecificInternalEnergy", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "HydroBase::press")
- call CCTK_VarIndex(dep_gfs(2), "HydroBase::eps")
- else
- call CCTK_VarIndex(indep_gfs(1), "HydroBase::rho")
- call CCTK_VarIndex(indep_gfs(2), "HydroBase::eps")
- call Util_TableSetString(ierr, table_handle, &
- "Pressure", &
- "Dependent variable names")
- call CCTK_VarIndex(dep_gfs(1), "HydroBase::press")
- end if
-
- call CCTK_FortranString(strlen, eosgeneral_indep_gfs, var_gfs)
- var_gfs = trim(var_gfs)
-
- strend = 1
-
- do var = 1, eosgeneral_n_indeps
-
- strend = scan(var_gfs, " ")
- if (strend > 0) then
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs(1:strend) )
- var_gfs = var_gfs(strend+1:)
- else
- if (var .ne. eosgeneral_n_indeps) then
- call CCTK_WARN(0, "Something wrong in the eos string")
- end if
- call CCTK_VarIndex( indep_gfs(n_base_indep_gfs+var), var_gfs )
- end if
-
- end do
-
- if (CCTK_EQUALS(GRHydro_eos_type, "Polytype")) then
- call Util_TableSetIntArray(ierr, table_handle, 1 + eosgeneral_n_indeps, &
- indep_gfs(1:1 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 2, dep_gfs, &
- "Dependent GFs")
- else
- call Util_TableSetIntArray(ierr, table_handle, 2 + eosgeneral_n_indeps, &
- indep_gfs(1:2 + eosgeneral_n_indeps), "Independent GFs")
- call Util_TableSetIntArray(ierr, table_handle, 1, dep_gfs, &
- "Dependent GFs")
- end if
-
- EOS_Prim2ConCellsCall = EOS_SetupCall(table_handle)
-
- end if
-
- deallocate(indep_gfs)
-
GRHydro_Startup = 0
end function GRHydro_Startup
diff --git a/src/make.code.defn b/src/make.code.defn
index 738ba0a..27c2ecb 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -15,7 +15,6 @@ SRCS = Utils.F90 \
GRHydro_EOS.c \
GRHydro_Flux.F90 \
GRHydro_FluxSplit.F90 \
- GRHydro_HLLE_EOSG.F90 \
GRHydro_HLLE.F90 \
GRHydro_Loop.F90 \
GRHydro_Minima.F90 \
@@ -49,7 +48,6 @@ SRCS = Utils.F90 \
GRHydro_Con2PrimM_polytype_pt.c \
GRHydro_EigenproblemM.F90 \
GRHydro_FluxM.F90 \
- GRHydro_HLLE_EOSGM.F90 \
GRHydro_HLLEM.F90 \
GRHydro_PPMM.F90 \
GRHydro_Prim2ConM.F90 \
@@ -59,9 +57,6 @@ SRCS = Utils.F90 \
GRHydro_UpdateMaskM.F90 \
GRHydro_UtilsM.F90
-### GRHydro_RegisterGZM.cc \
-### GRHydro_RegisterVarsM.cc \
-### GRHydro_Weights.c \
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-01 22:03:49 +0000