[Commits] [svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/src/ (Rev. 169)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Mon Nov 1 01:43:05 CDT 2010
User: cott
Date: 2010/11/01 01:43 AM
Modified:
/branches/hot_and_MHD_temp_dev/src/
GRHydro_Con2Prim.F90, GRHydro_Eigenproblem.F90, GRHydro_HLLEPoly.F90, GRHydro_ParamCheck.F90, GRHydro_Prim2Con.F90, GRHydro_RiemannSolve.F90, GRHydro_UpdateMask.F90
Log:
* bunch of small bug fixes in the 'hot' part of the code
File Changes:
Directory: /branches/hot_and_MHD_temp_dev/src/
==============================================
File [modified]: GRHydro_Con2Prim.F90
Delta lines: +2 -3
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_Con2Prim.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_Con2Prim.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -580,9 +580,9 @@
#if USE_EOS_OMNI
call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
rho,epsilon,temp,ye,xpress,keyerr,anyerr)
- pold = max(1.d-10,xpress)
+ pold = max(1.d-15,xpress)
#else
- pold = max(1.d-10,EOS_Pressure(handle, rho, epsilon))
+ pold = max(1.d-15,EOS_Pressure(handle, rho, epsilon))
#endif
@@ -615,7 +615,6 @@
#if USE_EOS_OMNI
call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
rho,epsilon,temp,ye,xpress,keyerr,anyerr)
-
f = pold - xpress
#else
f = pold - EOS_Pressure(handle, rho, epsilon)
File [modified]: GRHydro_Eigenproblem.F90
Delta lines: +73 -0
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_Eigenproblem.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_Eigenproblem.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -136,6 +136,79 @@
end subroutine eigenvalues
+
+subroutine eigenvalues_hot(handle,rho,velx,vely,velz,eps, &
+ temp,ye,w_lorentz,lam,gxx,gxy,gxz,gyy,gyz,gzz,u,alp,beta)
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL rho,velx,vely,velz,eps,w_lorentz
+ CCTK_REAL lam(5)
+ CCTK_REAL gxx,gxy,gxz,gyy,gyz,gzz
+ CCTK_REAL alp,beta,u
+ CCTK_REAL temp,ye
+
+ CCTK_REAL cs2,one,two
+ CCTK_REAL vlowx,vlowy,vlowz,v2,w
+ CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta
+ CCTK_INT handle
+ CCTK_REAL dpdrho,dpdeps,press
+
+#if !USE_EOS_OMNI
+#ifdef _EOS_BASE_INC_
+#undef _EOS_BASE_INC_
+#endif
+#include "EOS_Base.inc"
+#endif
+
+#if USE_EOS_OMNI
+! begin EOS Omni vars
+ integer :: n = 1
+ integer :: keytemp = 0
+ integer :: anyerr = 0
+ integer :: keyerr(1) = 0
+ real*8 :: xpress = 0.0d0
+ real*8 :: xeps = 0.0d0
+! end EOS Omni vars
+#endif
+
+ one = 1.0d0
+ two = 2.0d0
+
+!!$ Set required fluid quantities
+ call EOS_Omni_cs2(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho,eps,temp,ye,cs2,keyerr,anyerr)
+
+ vlowx = gxx*velx + gxy*vely + gxz*velz
+ vlowy = gxy*velx + gyy*vely + gyz*velz
+ vlowz = gxz*velx + gyz*vely + gzz*velz
+ v2 = vlowx*velx + vlowy*vely + vlowz*velz
+
+ w = w_lorentz
+
+!!$ Calculate eigenvalues
+
+ lam1 = velx - beta/alp
+ lam2 = velx - beta/alp
+ lam3 = velx - beta/alp
+ lamp_nobeta = (velx*(one-cs2) + sqrt(cs2*(one-v2)*&
+ (u*(one-v2*cs2) - velx**2*(one-cs2))))/(one-v2*cs2)
+ lamm_nobeta = (velx*(one-cs2) - sqrt(cs2*(one-v2)*&
+ (u*(one-v2*cs2) - velx**2*(one-cs2))))/(one-v2*cs2)
+
+ lamp = lamp_nobeta - beta/alp
+ lamm = lamm_nobeta - beta/alp
+
+ lam(1) = lamm
+ lam(2) = lam1
+ lam(3) = lam2
+ lam(4) = lam3
+ lam(5) = lamp
+
+
+end subroutine eigenvalues_hot
+
/*@@
@routine eigenproblem
@date Sat Jan 26 01:27:59 2002
File [modified]: GRHydro_HLLEPoly.F90
Delta lines: +123 -48
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_HLLEPoly.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_HLLEPoly.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -46,6 +46,7 @@
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
+ CCTK_REAL :: xtemp
CCTK_INT :: type_bits, trivial, not_trivial
@@ -179,22 +180,47 @@
!!$ Eigenvalues and fluxes either side of the cell interface
if (flux_direction == 1) then
- call eigenvalues(GRHydro_eos_handle,&
- rhominus(i+xoffset,j+yoffset,k+zoffset),&
- velxminus(i+xoffset,j+yoffset,k+zoffset),&
- velyminus(i+xoffset,j+yoffset,k+zoffset),&
- velzminus(i+xoffset,j+yoffset,k+zoffset),&
- epsminus(i+xoffset,j+yoffset,k+zoffset),&
- w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
- lamminus,gxxh,gxyh,gxzh,gyyh,&
- gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
- velxplus(i,j,k),velyplus(i,j,k),&
- velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
- lamplus,gxxh,gxyh,gxzh,gyyh,&
- gyzh,gzzh,&
- usendh,avg_alp,avg_beta)
+ if(evolve_temp.ne.1) then
+ call eigenvalues(GRHydro_eos_handle,&
+ rhominus(i+xoffset,j+yoffset,k+zoffset),&
+ velxminus(i+xoffset,j+yoffset,k+zoffset),&
+ velyminus(i+xoffset,j+yoffset,k+zoffset),&
+ velzminus(i+xoffset,j+yoffset,k+zoffset),&
+ epsminus(i+xoffset,j+yoffset,k+zoffset),&
+ w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gxxh,gxyh,gxzh,gyyh,&
+ gyzh,gzzh,&
+ usendh,avg_alp,avg_beta)
+ call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
+ velxplus(i,j,k),velyplus(i,j,k),&
+ velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
+ lamplus,gxxh,gxyh,gxzh,gyyh,&
+ gyzh,gzzh,&
+ usendh,avg_alp,avg_beta)
+ else
+ xtemp = temperature(i,j,k)
+ call eigenvalues_hot(GRHydro_eos_handle,&
+ rhominus(i+xoffset,j+yoffset,k+zoffset),&
+ velxminus(i+xoffset,j+yoffset,k+zoffset),&
+ velyminus(i+xoffset,j+yoffset,k+zoffset),&
+ velzminus(i+xoffset,j+yoffset,k+zoffset),&
+ epsminus(i+xoffset,j+yoffset,k+zoffset),&
+ xtemp,&
+ y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
+ w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gxxh,gxyh,gxzh,gyyh,&
+ gyzh,gzzh,&
+ usendh,avg_alp,avg_beta)
+ xtemp = temperature(i,j,k)
+ call eigenvalues_hot(GRHydro_eos_handle,rhoplus(i,j,k),&
+ velxplus(i,j,k),velyplus(i,j,k),&
+ velzplus(i,j,k),epsplus(i,j,k), &
+ xtemp,y_e_plus(i,j,k),&
+ w_lorentzplus(i,j,k),&
+ lamplus,gxxh,gxyh,gxzh,gyyh,&
+ gyzh,gzzh,&
+ usendh,avg_alp,avg_beta)
+ endif
call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),&
fplus(1),fplus(2),fplus(3),fplus(4),&
fplus(5),velxplus(i,j,k),pressplus(i,j,k),&
@@ -206,22 +232,47 @@
pressminus(i+xoffset,j+yoffset,k+zoffset),&
avg_det,avg_alp,avg_beta)
else if (flux_direction == 2) then
- call eigenvalues(GRHydro_eos_handle,&
- rhominus(i+xoffset,j+yoffset,k+zoffset),&
- velyminus(i+xoffset,j+yoffset,k+zoffset),&
- velzminus(i+xoffset,j+yoffset,k+zoffset),&
- velxminus(i+xoffset,j+yoffset,k+zoffset),&
- epsminus(i+xoffset,j+yoffset,k+zoffset),&
- w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
- lamminus,gyyh,gyzh,gxyh,gzzh,&
- gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
- velyplus(i,j,k),velzplus(i,j,k),&
- velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
- lamplus,gyyh,gyzh,gxyh,gzzh,&
- gxzh,gxxh,&
- usendh,avg_alp,avg_beta)
+ if(evolve_temp.ne.1) then
+ call eigenvalues(GRHydro_eos_handle,&
+ rhominus(i+xoffset,j+yoffset,k+zoffset),&
+ velyminus(i+xoffset,j+yoffset,k+zoffset),&
+ velzminus(i+xoffset,j+yoffset,k+zoffset),&
+ velxminus(i+xoffset,j+yoffset,k+zoffset),&
+ epsminus(i+xoffset,j+yoffset,k+zoffset),&
+ w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gyyh,gyzh,gxyh,gzzh,&
+ gxzh,gxxh,&
+ usendh,avg_alp,avg_beta)
+ call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
+ velyplus(i,j,k),velzplus(i,j,k),&
+ velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
+ lamplus,gyyh,gyzh,gxyh,gzzh,&
+ gxzh,gxxh,&
+ usendh,avg_alp,avg_beta)
+ else
+ xtemp = temperature(i,j,k)
+ call eigenvalues_hot(GRHydro_eos_handle,&
+ rhominus(i+xoffset,j+yoffset,k+zoffset),&
+ velyminus(i+xoffset,j+yoffset,k+zoffset),&
+ velzminus(i+xoffset,j+yoffset,k+zoffset),&
+ velxminus(i+xoffset,j+yoffset,k+zoffset),&
+ epsminus(i+xoffset,j+yoffset,k+zoffset),&
+ xtemp,&
+ y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
+ w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gyyh,gyzh,gxyh,gzzh,&
+ gxzh,gxxh,&
+ usendh,avg_alp,avg_beta)
+ xtemp = temperature(i,j,k)
+ call eigenvalues_hot(GRHydro_eos_handle,rhoplus(i,j,k),&
+ velyplus(i,j,k),velzplus(i,j,k),&
+ velxplus(i,j,k),epsplus(i,j,k),&
+ xtemp,y_e_plus(i,j,k),&
+ w_lorentzplus(i,j,k),&
+ lamplus,gyyh,gyzh,gxyh,gzzh,&
+ gxzh,gxxh,&
+ usendh,avg_alp,avg_beta)
+ endif
call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),&
fplus(1),fplus(3),fplus(4),fplus(2),&
fplus(5),velyplus(i,j,k),pressplus(i,j,k),&
@@ -233,22 +284,46 @@
pressminus(i+xoffset,j+yoffset,k+zoffset),&
avg_det,avg_alp,avg_beta)
else if (flux_direction == 3) then
- call eigenvalues(GRHydro_eos_handle,&
- rhominus(i+xoffset,j+yoffset,k+zoffset),&
- velzminus(i+xoffset,j+yoffset,k+zoffset),&
- velxminus(i+xoffset,j+yoffset,k+zoffset),&
- velyminus(i+xoffset,j+yoffset,k+zoffset),&
- epsminus(i+xoffset,j+yoffset,k+zoffset),&
- w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
- lamminus,gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
- call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
- velzplus(i,j,k),velxplus(i,j,k),&
- velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
- lamplus,gzzh,gxzh,gyzh,&
- gxxh,gxyh,gyyh,&
- usendh,avg_alp,avg_beta)
+ if(evolve_temp.ne.1) then
+ call eigenvalues(GRHydro_eos_handle,&
+ rhominus(i+xoffset,j+yoffset,k+zoffset),&
+ velzminus(i+xoffset,j+yoffset,k+zoffset),&
+ velxminus(i+xoffset,j+yoffset,k+zoffset),&
+ velyminus(i+xoffset,j+yoffset,k+zoffset),&
+ epsminus(i+xoffset,j+yoffset,k+zoffset),&
+ w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gzzh,gxzh,gyzh,&
+ gxxh,gxyh,gyyh,&
+ usendh,avg_alp,avg_beta)
+ call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),&
+ velzplus(i,j,k),velxplus(i,j,k),&
+ velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),&
+ lamplus,gzzh,gxzh,gyzh,&
+ gxxh,gxyh,gyyh,&
+ usendh,avg_alp,avg_beta)
+ else
+ xtemp = temperature(i,j,k)
+ call eigenvalues_hot(GRHydro_eos_handle,&
+ rhominus(i+xoffset,j+yoffset,k+zoffset),&
+ velzminus(i+xoffset,j+yoffset,k+zoffset),&
+ velxminus(i+xoffset,j+yoffset,k+zoffset),&
+ velyminus(i+xoffset,j+yoffset,k+zoffset),&
+ epsminus(i+xoffset,j+yoffset,k+zoffset),&
+ xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),&
+ w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),&
+ lamminus,gzzh,gxzh,gyzh,&
+ gxxh,gxyh,gyyh,&
+ usendh,avg_alp,avg_beta)
+ xtemp = temperature(i,j,k)
+ call eigenvalues_hot(GRHydro_eos_handle,rhoplus(i,j,k),&
+ velzplus(i,j,k),velxplus(i,j,k),&
+ velyplus(i,j,k),epsplus(i,j,k),&
+ xtemp,y_e_plus(i,j,k),&
+ w_lorentzplus(i,j,k),&
+ lamplus,gzzh,gxzh,gyzh,&
+ gxxh,gxyh,gyyh,&
+ usendh,avg_alp,avg_beta)
+ endif
call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),&
fplus(1),fplus(4),fplus(2),fplus(3),&
fplus(5),velzplus(i,j,k),pressplus(i,j,k),avg_det,&
File [modified]: GRHydro_ParamCheck.F90
Delta lines: +1 -1
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_ParamCheck.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_ParamCheck.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -98,7 +98,7 @@
evolve_Y_e = 0
endif
- if (CCTK_EQUALS(temperature_and_entropy_evolution_method,"GRHydro")) then
+ if (CCTK_EQUALS(temp_and_ent_evolution_method,"GRHydro")) then
evolve_temp = 1
else
evolve_temp = 0
File [modified]: GRHydro_Prim2Con.F90
Delta lines: +0 -3
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -141,9 +141,6 @@
endif
-
-
-
end subroutine primitive2conservative
/*@@
File [modified]: GRHydro_RiemannSolve.F90
Delta lines: +1 -1
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_RiemannSolve.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_RiemannSolve.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -39,7 +39,7 @@
CCTK_INT :: i,j,k
if (CCTK_EQUALS(riemann_solver,"HLLE")) then
-
+
call GRHydro_HLLE(CCTK_PASS_FTOF)
if (evolve_tracer .ne. 0) then
File [modified]: GRHydro_UpdateMask.F90
Delta lines: +8 -8
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_UpdateMask.F90 2010-11-01 05:49:09 UTC (rev 168)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_UpdateMask.F90 2010-11-01 06:43:05 UTC (rev 169)
@@ -297,10 +297,10 @@
vely(i,j,k) = 0.0d0
velz(i,j,k) = 0.0d0
#if USE_EOS_OMNI
- call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
- rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
- call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
- rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
+ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr)
+ call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr)
#else
press(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps(i,j,k))
eps(i,j,k) = EOS_SpecificIntEnergy(eos_handle, rho_min, press(i,j,k))
@@ -321,9 +321,9 @@
vely_p(i,j,k) = 0.0d0
velz_p(i,j,k) = 0.0d0
#if USE_EOS_OMNI
- call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr)
- call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr)
#else
press_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p(i,j,k))
@@ -346,9 +346,9 @@
vely_p_p(i,j,k) = 0.0d0
velz_p_p(i,j,k) = 0.0d0
#if USE_EOS_OMNI
- call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr)
- call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr)
#else
press_p_p(i,j,k) = EOS_Pressure(eos_handle, rho_min, eps_p_p(i,j,k))
More information about the Commits
mailing list