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))