[Commits] [svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/ (Rev. 176)

User: cott Date: 2010/11/16 09:19 AM

Modified: /branches/hot_and_MHD_temp_dev/ schedule.ccl /branches/hot_and_MHD_temp_dev/src/ GRHydro_Boundaries.F90, GRHydro_Con2Prim.F90, GRHydro_Prim2Con.F90, GRHydro_Source.F90, GRHydro_UpdateMask.F90, Utils.F90

Log: * a number of small modifications having to do with boundary/atmosphere treatment when the hot EOS is being used.

File Changes:

Directory: /branches/hot_and_MHD_temp_dev/ ==========================================

File [modified]: schedule.ccl Delta lines: +11 -0 =================================================================== --- branches/hot_and_MHD_temp_dev/schedule.ccl 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/schedule.ccl 2010-11-16 15:19:43 UTC (rev 176) @@ -877,6 +877,17 @@ LANG: Fortran } "Calculate current refinement level"

+#debug +### schedule GRHydro_Debug IN MoL_PostRHS AFTER dissipation_add +### { +### LANG: Fortran +### } "debug" +### +### schedule GRHydro_Debug AS guender IN MoL_Step AFTER MoL_Add +### { +### LANG: Fortran +### } "debug" + schedule GRHydro_RefinementLevel AT CCTK_PostStep { LANG: Fortran

Directory: /branches/hot_and_MHD_temp_dev/src/ ==============================================

File [modified]: GRHydro_Boundaries.F90 Delta lines: +2 -2 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Boundaries.F90 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Boundaries.F90 2010-11-16 15:19:43 UTC (rev 176) @@ -66,12 +66,12 @@ call SetCartSymVN(ierr, cctkGH, sym, "GRHydro::GRHydro_C2P_failed")

!!$ handle multiple tracer variables - if(evolve_tracer.eq.1) then + if(evolve_tracer.ne.0) then call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_tracers") call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::GRHydro_cons_tracers") endif

- if(CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then + if(evolve_y_e.ne.0) then call SetCartSymGN(ierr, cctkGH, sym, "HydroBase::Y_e") call SetCartSymGN(ierr, cctkGH, sym, "GRHydro::Y_e_con") endif

File [modified]: GRHydro_Con2Prim.F90 Delta lines: +18 -1 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Con2Prim.F90 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Con2Prim.F90 2010-11-16 15:19:43 UTC (rev 176) @@ -230,7 +230,6 @@ end do !$OMP END PARALLEL DO

- return

end subroutine Conservative2Primitive @@ -502,6 +501,10 @@ real*8 :: xpress = 0.0d0 ! end EOS Omni vars

+! set pmin and epsmin to something sensible: + pmin = 1.0d-28 + epsmin = 1.0e-5 + !!$ Undensitize the variables

udens = dens /sqrt(det) @@ -516,7 +519,21 @@ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) pold = max(1.d-15,xpress) + ! error handling + if(anyerr.ne.0) then + call CCTK_WARN(1,"EOS error in c2p 0") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif

+ !!$ Check that the variables have a chance of being physical

if( (utau + pold + udens)**2 - s2 .le. 0.0d0) then

File [modified]: GRHydro_Prim2Con.F90 Delta lines: +25 -6 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-16 15:19:43 UTC (rev 176) @@ -113,7 +113,8 @@ ! eps is reconstructed. Hence, we do not update the ! variable 'temperature' in prim2con at the interfaces xtemp = temperature(i,j,k) - call prim2con_hot(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,& + call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxl,gxyl,gxzl,gyyl,& gyzl,gzzl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& @@ -126,7 +127,8 @@ ! eps is reconstructed. Hence, we do not update the ! variable 'temperature' in prim2con at the interfaces xtemp = temperature(i,j,k) - call prim2con_hot(GRHydro_eos_handle, gxxr,gxyr,gxzr,& + call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxr,gxyr,gxzr,& gyyr,gyzr,gzzr, & avg_detr, densplus(i,j,k),sxplus(i,j,k),& syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& @@ -201,7 +203,8 @@ end subroutine prim2con

-subroutine prim2con_hot(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & +subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & + x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, & temp,ye)

@@ -212,8 +215,9 @@ CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& deps, dpress, w, vlowx, vlowy, vlowz - CCTK_REAL :: temp,ye - CCTK_INT :: handle + CCTK_REAL :: temp,ye, x, y, z + CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk + character(len=512) warnline

! begin EOS Omni vars integer :: n = 1 @@ -231,6 +235,19 @@

call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) + ! error handling + if(anyerr.ne.0) then + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif

vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz @@ -271,6 +288,7 @@

CCTK_INT :: i, j, k CCTK_REAL :: det + character(len=512) :: warnline

if(evolve_Y_e.ne.0) then !$OMP PARALLEL DO PRIVATE(i, j) @@ -312,7 +330,8 @@ det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))

- call prim2con_hot(GRHydro_eos_handle,gxx(i,j,k),& + call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,i,j,k,& + x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),& gxy(i,j,k),gxz(i,j,k),& gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&

File [modified]: GRHydro_Source.F90 Delta lines: +1 -1 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Source.F90 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Source.F90 2010-11-16 15:19:43 UTC (rev 176) @@ -101,7 +101,7 @@ end if

if (evolve_Y_e .ne. 0) then - Y_e_con_rhs = 0.0d0 + y_e_con_rhs = 0.0d0 endif

!!$ Set up the array for checking the order. We switch to second order

File [modified]: GRHydro_UpdateMask.F90 Delta lines: +89 -36 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_UpdateMask.F90 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_UpdateMask.F90 2010-11-16 15:19:43 UTC (rev 176) @@ -123,6 +123,10 @@ cons_tracerrhs = 0.0d0 endif

+ if (evolve_y_e .ne. 0) then + y_e_con_rhs = 0.0d0 + endif + atmosphere_mask = 0

call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere") @@ -308,23 +312,40 @@

if (rho(i,j,k) .le. rho_min) then rho(i,j,k) = rho_min + det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \ + gyy(i,j,k), gyz(i,j,k), gzz(i,j,k)) + velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.0d0

- 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) + if(evolve_temper.ne.0) then + ! set the temperature to be relatively low + temperature(i,j,k) = 0.1d0 + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr)

- det = SPATIAL_DETERMINANT(gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), \ - gyy(i,j,k), gyz(i,j,k), gzz(i,j,k)) - call prim2conpolytype(eos_handle, & - gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & - gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & - dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & - tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & - velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) + call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & + det,dens(i,j,k),scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & + velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k),& + temperature(i,j,k),y_e(i,j,k)) + else + 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) + call prim2conpolytype(eos_handle, & + gxx(i,j,k), gxy(i,j,k), gxz(i,j,k), & + gyy(i,j,k), gyz(i,j,k), gzz(i,j,k), det, & + dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + tau(i,j,k), rho(i,j,k), velx(i,j,k), vely(i,j,k), & + velz(i,j,k), eps(i,j,k), press(i,j,k), w_lorentz(i,j,k)) + endif end if if (timelevels .gt. 1) then if (rho_p(i,j,k) .le. rho_min) then @@ -333,19 +354,37 @@ vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0

- 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(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) - det = SPATIAL_DETERMINANT(gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), \ gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k)) - call prim2conpolytype(eos_handle, & - gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), & - gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, & - dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & - tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & - velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k)) + + if(evolve_temper.ne.0) then + ! set the temperature to be relatively low + temperature_p(i,j,k) = 0.1d0 + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),eps_p(i,j,k),temperature_p(i,j,k),y_e_p(i,j,k),& + press_p(i,j,k),keyerr,anyerr) + + call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p(i,j,k),gxy_p(i,j,k),& + gxz_p(i,j,k),gyy_p(i,j,k),gyz_p(i,j,k),gzz_p(i,j,k), & + det,dens_p(i,j,k),scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & + tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & + velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k),& + temperature_p(i,j,k),y_e_p(i,j,k)) + else + 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(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) + call prim2conpolytype(eos_handle, & + gxx_p(i,j,k), gxy_p(i,j,k), gxz_p(i,j,k), & + gyy_p(i,j,k), gyz_p(i,j,k), gzz_p(i,j,k), det, & + dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & + tau_p(i,j,k), rho_p(i,j,k), velx_p(i,j,k), vely_p(i,j,k), & + velz_p(i,j,k), eps_p(i,j,k), press_p(i,j,k), w_lorentz_p(i,j,k)) + endif + endif end if if (timelevels .gt. 2) then @@ -354,20 +393,35 @@ velx_p_p(i,j,k) = 0.0d0 vely_p_p(i,j,k) = 0.0d0 velz_p_p(i,j,k) = 0.0d0 - - 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(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) - det = SPATIAL_DETERMINANT(gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), \ gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k)) - call prim2conpolytype(eos_handle, & - gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), & - gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, & - dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & - tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & - velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k)) + + if(evolve_temper.ne.0) then + ! set the temperature to be relatively low + temperature_p_p(i,j,k) = 0.1d0 + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),& + press_p_p(i,j,k),keyerr,anyerr) + call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxx_p_p(i,j,k),gxy_p_p(i,j,k),& + gxz_p_p(i,j,k),gyy_p_p(i,j,k),gyz_p_p(i,j,k),gzz_p_p(i,j,k), & + det,dens_p_p(i,j,k),scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & + tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & + velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k),& + temperature_p_p(i,j,k),y_e_p_p(i,j,k)) + else + 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(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) + call prim2conpolytype(eos_handle, & + gxx_p_p(i,j,k), gxy_p_p(i,j,k), gxz_p_p(i,j,k), & + gyy_p_p(i,j,k), gyz_p_p(i,j,k), gzz_p_p(i,j,k), det, & + dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & + tau_p_p(i,j,k), rho_p_p(i,j,k), velx_p_p(i,j,k), vely_p_p(i,j,k), & + velz_p_p(i,j,k), eps_p_p(i,j,k), press_p_p(i,j,k), w_lorentz_p_p(i,j,k)) + endif endif endif

@@ -375,6 +429,5 @@ end do end do

- end subroutine GRHydro_InitialAtmosphereReset

File [modified]: Utils.F90 Delta lines: +31 -0 =================================================================== --- branches/hot_and_MHD_temp_dev/src/Utils.F90 2010-11-14 05:28:48 UTC (rev 175) +++ branches/hot_and_MHD_temp_dev/src/Utils.F90 2010-11-16 15:19:43 UTC (rev 176) @@ -27,6 +27,37 @@ @endhistory @@*/

+subroutine GRHydro_Debug(CCTK_ARGUMENTS) + + implicit none + DECLARE_CCTK_ARGUMENTS + integer loc(3) + + loc = minloc(y_e_con/dens) +#if 0 + write(6,*) minval(y_e_con/dens) + write(6,*) y_e_con_rhs(loc(1),loc(2),loc(3)),densrhs(loc(1),loc(2),loc(3)) + write(6,*) y_e_con(loc(1),loc(2),loc(3)),dens(loc(1),loc(2),loc(3)) + write(6,*) loc + write(6,*) dens(4,5,5),y_e(4,5,5),y_e_con(4,5,5) + write(6,*) x(loc(1),loc(2),loc(3)),y(loc(1),loc(2),loc(3)),z(loc(1),loc(2),loc(3)) +#endif + loc(1) = 4 + loc(2) = 5 + loc(3) = 4 + write(6,*) minval(y_e_con/dens) + write(6,*) y_e_con_rhs(loc(1),loc(2),loc(3)),densrhs(loc(1),loc(2),loc(3)) + write(6,*) y_e_con(loc(1),loc(2),loc(3)),dens(loc(1),loc(2),loc(3)) + write(6,*) loc + write(6,*) dens(4,5,5),y_e(4,5,5),y_e_con(4,5,5) + write(6,*) x(loc(1),loc(2),loc(3)),y(loc(1),loc(2),loc(3)),z(loc(1),loc(2),loc(3)) + + + call flush(6) + +end subroutine GRHydro_Debug + + subroutine GRHydro_RefinementLevel(CCTK_ARGUMENTS)

implicit none