[Commits] [svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/src/ (Rev. 167)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Mon Nov 1 00:33:27 CDT 2010
User: cott
Date: 2010/11/01 12:33 AM
Modified:
/branches/hot_and_MHD_temp_dev/src/
GRHydro_Prim2Con.F90
Log:
* add hot p2c routines
File Changes:
Directory: /branches/hot_and_MHD_temp_dev/src/
==============================================
File [modified]: GRHydro_Prim2Con.F90
Delta lines: +188 -53
===================================================================
--- branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-10-31 16:19:38 UTC (rev 166)
+++ branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-01 05:33:26 UTC (rev 167)
@@ -46,44 +46,104 @@
integer :: i, j, k
CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
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))
+ CCTK_REAL :: xtemp
- avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)
- avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)
+ if(evolve_temp.ne.1) then
+ 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 prim2con(GRHydro_eos_handle, 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 prim2con(GRHydro_eos_handle, 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
+ else
+ 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 prim2con(GRHydro_eos_handle, 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 prim2con(GRHydro_eos_handle, 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))
+ ! we do not have plus/minus vars for temperature since
+ ! 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,&
+ 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),&
+ xtemp,y_e_minus(i,j,k))
- end do
- end do
- end do
+ ! we do not have plus/minus vars for temperature since
+ ! 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,&
+ 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),xtemp, &
+ y_e_plus(i,j,k))
+
+ end do
+ end do
+ end do
+ endif
+
+
+
+
end subroutine primitive2conservative
/*@@
@@ -154,6 +214,59 @@
end subroutine prim2con
+subroutine prim2con_hot(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, &
+ dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, &
+ temp,ye)
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ 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
+
+#if !USE_EOS_OMNI
+#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
+
+ w = 1.d0 / sqrt(1.d0 - (gxx*dvelx*dvelx + gyy*dvely*dvely + gzz &
+ *dvelz*dvelz + 2*gxy*dvelx*dvely + 2*gxz*dvelx *dvelz + 2*gyz&
+ *dvely*dvelz))
+
+#if USE_EOS_OMNI
+ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&
+ drho,deps,temp,ye,dpress,keyerr,anyerr)
+#else
+ dpress = EOS_Pressure(handle, drho, deps)
+#endif
+
+ vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz
+ vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
+ vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
+
+ ddens = sqrt(det) * drho * w
+ dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx
+ dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy
+ dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz
+ dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens
+
+end subroutine prim2con_hot
+
+
/*@@
@routine Primitive2ConservativeCells
@date Sun Mar 10 21:16:20 2002
@@ -180,26 +293,7 @@
CCTK_INT :: i, j, k
CCTK_REAL :: det
-
- !$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
-
- 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(GRHydro_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),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
- tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
- eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
- end do
- end do
- end do
- !$OMP END PARALLEL DO
-
if(evolve_Y_e.ne.0) then
!$OMP PARALLEL DO PRIVATE(i, j)
do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1
@@ -212,7 +306,48 @@
!$OMP END PARALLEL DO
endif
+ if(evolve_temp.ne.1) 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
+
+ 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(GRHydro_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),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
+ tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),&
+ eps(i,j,k),press(i,j,k),w_lorentz(i,j,k))
+ end do
+ end do
+ end do
+ !$OMP END PARALLEL DO
+ else
+ !$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
+
+ 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),&
+ 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),&
+ 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))
+ end do
+ end do
+ end do
+ !$OMP END PARALLEL DO
+ endif
+
+
end subroutine Primitive2ConservativeCells
More information about the Commits
mailing list