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