diff options
Diffstat (limited to 'src/GRHydro_Prim2ConAM.F90')
| -rw-r--r-- | src/GRHydro_Prim2ConAM.F90 | 114 |
1 files changed, 57 insertions, 57 deletions
diff --git a/src/GRHydro_Prim2ConAM.F90 b/src/GRHydro_Prim2ConAM.F90 index 7bda654..2aee515 100644 --- a/src/GRHydro_Prim2ConAM.F90 +++ b/src/GRHydro_Prim2ConAM.F90 @@ -54,8 +54,8 @@ subroutine primitive2conservativeAM(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 + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_sdetl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_sdetr !CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,& ! g11r,g12r,g13r,g22r,g23r,g33r CCTK_REAL :: xtemp(1) @@ -76,7 +76,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& + !$OMP PARALLEL DO PRIVATE(k,j,i,avg_sdetl,avg_sdetr,& !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,& !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -96,11 +96,11 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) 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) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) call prim2conAM(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,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),& @@ -108,7 +108,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k)) call prim2conAM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, 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),& @@ -122,7 +122,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) else - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -142,8 +142,8 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) 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) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr)) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -161,7 +161,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) !$OMP END CRITICAL call prim2conAM_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),& + avg_sdetl,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),& @@ -177,7 +177,7 @@ subroutine primitive2conservativeAM(CCTK_ARGUMENTS) call prim2conAM_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),& + avg_sdetr, 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),& @@ -208,15 +208,17 @@ end subroutine primitive2conservativeAM @@*/ -subroutine prim2conAM_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, dBvcx, dBvcy, dBvcz, w, temp, ye) +subroutine prim2conAM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, & + gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & + dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, & + deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz,& deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz @@ -333,23 +335,23 @@ subroutine prim2conAM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gx 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 - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((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 + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens !!$OMP CRITICAL - !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w + !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sdet, w:)', ddens,sdet,w !call CCTK_WARN(1, NaN_WarnLine) !!$OMP END CRITICAL end subroutine prim2conAM_hot -subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & +subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none @@ -357,7 +359,7 @@ subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz,& deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz @@ -403,14 +405,14 @@ subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & 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 - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((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 + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens end subroutine prim2conAM @@ -439,25 +441,24 @@ subroutine Primitive2ConservativeCellsAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: maxtau0 maxtau0 = -1.0d60 if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(k,j,i,det), REDUCTION(MAX:maxtau0) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet), REDUCTION(MAX:maxtau0) 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)) + sdet = sdetg(i,j,k) call prim2conAM(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),& + sdet, 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),Bvecx(i,j,k), & @@ -483,19 +484,18 @@ subroutine Primitive2ConservativeCellsAM(CCTK_ARGUMENTS) !GRHydro_tau_min = GRHydro_tau_min * maxtau0 else - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) 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)) + sdet = sdetg(i,j,k) call prim2conAM_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),& + sdet, 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),Bvecx(i,j,k), & @@ -548,13 +548,13 @@ subroutine Prim2ConservativePolytypeAM(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 + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_sdetl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_sdetr ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces - !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_sdetl,& + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_sdetr) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) @@ -572,12 +572,12 @@ subroutine Prim2ConservativePolytypeAM(CCTK_ARGUMENTS) 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) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr)) call prim2conpolytypeAM(GRHydro_eos_handle, gxxl,gxyl,gxzl,& gyyl,gyzl,gzzl, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,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),& @@ -586,7 +586,7 @@ subroutine Prim2ConservativePolytypeAM(CCTK_ARGUMENTS) call prim2conpolytypeAM(GRHydro_eos_handle, gxxr,gxyr,gxzr,& gyyr,gyzr,gzzr, & - avg_detr,densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr,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),& @@ -620,14 +620,14 @@ end subroutine Prim2ConservativePolytypeAM @@*/ subroutine prim2conpolytypeAM(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, dsx, dsy, dsz, dtau, & + gzz, sdet, ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdet CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & w_tmp, w, vlowx, vlowy, vlowz @@ -680,14 +680,14 @@ subroutine prim2conpolytypeAM(handle, gxx, gxy, gxz, gyy, gyz, & 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 - & + ddens = sdet * drho * w + dsx = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & ab0*blowx) - dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + dsy = sdet * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & ab0*blowy) - dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + dsz = sdet * ((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 + dtau = sdet * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens end subroutine prim2conpolytypeAM @@ -717,19 +717,19 @@ subroutine Primitive2ConservativePolyCellsAM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) 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)) + sdet = sdetg(i,j,k) + call prim2conpolytypeAM(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),& + sdet, 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),Bvecx(i,j,k), Bvecy(i,j,k), & |