/*@@ @file primitive2conservative @date Aug 31, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc Primitive to conservative routine @enddesc @@*/ #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "cctk_Functions.h" #include "GRHydro_Macros.h" #define velx(i,j,k) vel(i,j,k,1) #define vely(i,j,k) vel(i,j,k,2) #define velz(i,j,k) vel(i,j,k,3) #define sx(i,j,k) scon(i,j,k,1) #define sy(i,j,k) scon(i,j,k,2) #define sz(i,j,k) scon(i,j,k,3) #define Bvecx(i,j,k) Bvec(i,j,k,1) #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) #define DOT(x1,y1,z1,x2,y2,z2) ( DOTP(gxx,gxy,gxz,gyy,gyz,gzz,x1,y1,z1,x2,y2,z2) ) #define DOT2(x1,y1,z1) ( DOTP2(gxx,gxy,gxz,gyy,gyz,gzz,x1,y1,z1) ) #define DOTPT(x1,y1,z1,x2,y2,z2) ( DOTP(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt,x1,y1,z1,x2,y2,z2) ) #define DOTPT2(x1,y1,z1) ( DOTP2(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt,x1,y1,z1) ) /*@@ @routine primitive2conservativeM @date Aug 31 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc Converts primitive to conserved variables for the boundary extended data. @enddesc @calls @calledby @history @endhistory @@*/ subroutine primitive2conservativeM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS 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)) avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) call prim2conM(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),& Bvecxminus(i,j,k),Bvecyminus(i,j,k),Bveczminus(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 prim2conM(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),& Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(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 end subroutine primitive2conservativeM /*@@ @routine prim2conM @date Aug 31, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc Converts from primitive to conservative at a single point @enddesc @calls @calledby @history @endhistory @@*/ subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & dsx, dsy, dsz, dtau , dBvcx, dBvcy, dBvcz, drho, dvelx, dvely, dvelz, deps, dpress, w) implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & drho, dvelx, dvely, dvelz,& deps, dpress, vlowx, vlowy, vlowz CCTK_REAL :: w CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz CCTK_INT :: handle character(len=256) NaN_WarnLine ! begin EOS Omni vars integer :: n, keytemp, anyerr, keyerr(1) real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 ! end EOS Omni vars w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1))) !!$ BEGIN: Check for NaN value if (w .ne. w) then !$OMP CRITICAL write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) !$OMP END CRITICAL endif !!$ END: Check for NaN value call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,xtemp,xye,dpress,keyerr,anyerr) vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) ab0=w*Bdotv b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & w*Bdotv*vlowx blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & w*Bdotv*vlowy blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz ddens = sqrt(det) * drho * w dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens end subroutine prim2conM /*@@ @routine Primitive2ConservativeCellsM @date Aug 31, 2010 @author @desc Wrapper function that converts primitive to conservative at the cell centres. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k CCTK_REAL :: det 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 prim2conM(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),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(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 if(evolve_Y_e.ne.0) 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 Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k) enddo enddo enddo !$OMP END PARALLEL DO endif end subroutine Primitive2ConservativeCellsM /*@@ @routine Prim2ConservativePolytypeM @date Aug 31, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke @desc Same as first routine, only for polytropes. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS integer :: i, j, k CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& !$OMP 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)) avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) call prim2conpolytypeM(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),& Bvecxminus(i,j,k),Bvecyminus(i,j,k),Bveczminus(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 prim2conpolytypeM(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),& Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(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)) if (densminus(i,j,k) .ne. densminus(i,j,k)) then !$OMP CRITICAL call CCTK_WARN(1, "NaN in densminus(i,j,k) (Prim2Con)") !$OMP END CRITICAL endif end do end do end do !$OMP END PARALLEL DO end subroutine Prim2ConservativePolytypeM /*@@ @routine prim2conpolytypeM @date Aug 31, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc Converts from primitive to conservative at a single point @enddesc @calls @calledby @history @endhistory @@*/ subroutine prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, & gzz, det, ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & drho, dvelx, dvely, dvelz, deps, dpress, w) implicit none DECLARE_CCTK_PARAMETERS CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, dBvcx, dBvcy, dBvcz, & drho, dvelx, dvely, dvelz,& deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet CCTK_INT :: handle CCTK_REAL :: Bdotv,ab0,b2,blowx,blowy,blowz character(len=256) NaN_WarnLine ! begin EOS Omni vars integer :: n, keytemp, anyerr, keyerr(1) real*8 :: xpress,xeps,xtemp,xye n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 ! end EOS Omni vars w_tmp = DOT2(dvelx,dvely,dvelz) if (w_tmp .ge. 1.d0) then ! In theory this should not happen, and even when accepting the fact ! that numerically it can, one might be tempted to set w to some large ! value in that case. However, this would lead to completely bogus ! and hard to trace wrong values below. There is no good value to ! choose in this case, but something small is probably the best of ! all bad choices. !$OMP CRITICAL write(NaN_WarnLine,'(a80,2g15.6)') 'Infinite Lorentz factor reset. rho, w_tmp: ', drho, w_tmp call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) !$OMP END CRITICAL w = 1.d-20 else w = 1.d0 / sqrt(1.d0 - w_tmp) endif call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,xeps,xtemp,xye,dpress,keyerr,anyerr) call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,xeps,xtemp,xye,dpress,deps,keyerr,anyerr) vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) ab0=w*Bdotv b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & w*Bdotv*vlowx blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & w*Bdotv*vlowy blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & w*Bdotv*vlowz ddens = sqrt(det) * drho * w dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & ab0*blowz) dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens end subroutine prim2conpolytypeM /*@@ @routine Primitive2ConservativePolyCellsM @date Aug 31, 2010 @author @desc Wrapper function that converts primitive to conservative at the cell centres. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k CCTK_REAL :: det 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 prim2conpolytypeM(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),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(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 if(evolve_Y_e.ne.0) 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 Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k) enddo enddo enddo !$OMP END PARALLEL DO endif end subroutine Primitive2ConservativePolyCellsM !!$ !!$ Prim2Con doesn't change for tracers with the addition of a B-field! !!$ /*@@ @routine primitive2conservativegeneralM @date Aug 31, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc Converts primitive to conserved variables for the boundary extended data. @enddesc @calls @calledby @history @endhistory @@*/ subroutine primitive2conservativegeneralM(CCTK_ARGUMENTS) USE GRHydro_Scalars implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, ierr CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl, & gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr, & vlowx, vlowy, vlowz ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallPlus) ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConBndCallMinus) !$OMP PARALLEL DO PRIVATE(i, j, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& !$OMP 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)) avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) call prim2conpolytypeM(-1, 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),& Bvecxminus(i,j,k),Bvecyminus(i,j,k),Bveczminus(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 prim2conpolytypeM(-1, 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),& Bvecxplus(i,j,k),Bvecyplus(i,j,k),Bveczplus(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)) if (densminus(i,j,k) .ne. densminus(i,j,k)) then !$OMP CRITICAL call CCTK_WARN(1, "NaN in densminus(i,j,k) (Prim2Con)") !$OMP END CRITICAL endif end do end do end do !$OMP END PARALLEL DO end subroutine primitive2conservativegeneralM /*@@ @routine Primitive2ConservativeCellsGeneralM @date Ag 31, 2010 @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke @desc Converts primitive to conserved variables for the boundary extended data. @enddesc @calls @calledby @history @endhistory @@*/ subroutine Primitive2ConservativeCellsGeneralM(CCTK_ARGUMENTS) USE GRHydro_Scalars implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, ierr CCTK_REAL :: det,gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt, & vlowx, vlowy, vlowz CCTK_REAL Bdotv,ab0,b2,blowx,blowy,blowz ierr = EOS_SetGFs(cctkGH, EOS_Prim2ConCellsCall) 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 gxxpt = gxx(i,j,k) gxypt = gxy(i,j,k) gxzpt = gxz(i,j,k) gyypt = gyy(i,j,k) gyzpt = gyz(i,j,k) gzzpt = gzz(i,j,k) det = SPATIAL_DETERMINANT(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt) w_lorentz(i,j,k)=1.d0/sqrt(1.d0-DOTPT2(velx(i,j,k),vely(i,j,k),velz(i,j,k))) vlowx = gxxpt * velx(i,j,k) + & gxypt * vely(i,j,k) + & gxzpt * velz(i,j,k) vlowy = gxypt * velx(i,j,k) + & gyypt * vely(i,j,k) + & gyzpt * velz(i,j,k) vlowz = gxzpt * velx(i,j,k) + & gyzpt * vely(i,j,k) + & gzzpt * velz(i,j,k) Bdotv=DOTPT(velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k)) ab0=w_lorentz(i,j,k)*Bdotv b2 = DOTPT2(Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k))/w_lorentz(i,j,k)**2+Bdotv**2 blowx = (gxx(i,j,k)*Bvecx(i,j,k) + gxy(i,j,k)*Bvecy(i,j,k) + gxz(i,j,k)*Bvecz(i,j,k))/ & w_lorentz(i,j,k) + w_lorentz(i,j,k)*Bdotv*vlowx blowy = (gxy(i,j,k)*Bvecx(i,j,k) + gyy(i,j,k)*Bvecy(i,j,k) + gyz(i,j,k)*Bvecz(i,j,k))/ & w_lorentz(i,j,k) + w_lorentz(i,j,k)*Bdotv*vlowy blowz = (gxz(i,j,k)*Bvecx(i,j,k) + gyz(i,j,k)*Bvecy(i,j,k) + gzz(i,j,k)*Bvecz(i,j,k))/ & w_lorentz(i,j,k) + w_lorentz(i,j,k)*Bdotv*vlowz dens(i,j,k) = sqrt(det) * rho(i,j,k) * w_lorentz(i,j,k) sx(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* & w_lorentz(i,j,k)*w_lorentz(i,j,k) * vlowx - ab0*blowx) sy(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* & w_lorentz(i,j,k)*w_lorentz(i,j,k) * vlowy - ab0*blowy) sz(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* & w_lorentz(i,j,k)*w_lorentz(i,j,k) * vlowz - ab0*blowz) tau(i,j,k) = sqrt(det) * ((rho(i,j,k)*(1+eps(i,j,k))+press(i,j,k)+b2)* & w_lorentz(i,j,k)*w_lorentz(i,j,k) - press(i,j,k)-b2/2.0-ab0**2) - dens(i,j,k) end do end do end do end subroutine Primitive2ConservativeCellsGeneralM