diff options
Diffstat (limited to 'src/GRHydro_Prim2ConM.F90')
| -rw-r--r-- | src/GRHydro_Prim2ConM.F90 | 138 |
1 files changed, 69 insertions, 69 deletions
diff --git a/src/GRHydro_Prim2ConM.F90 b/src/GRHydro_Prim2ConM.F90 index cb459ee..a473163 100644 --- a/src/GRHydro_Prim2ConM.F90 +++ b/src/GRHydro_Prim2ConM.F90 @@ -61,8 +61,8 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr CCTK_REAL :: xtemp(1) character(len=256) NaN_WarnLine ! save memory when MP is not used @@ -96,7 +96,7 @@ subroutine primitive2conservativeM(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 g11l,g12l,g13l,g22l,g23l,g33l,& !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -116,11 +116,11 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2conM(GRHydro_eos_handle, g11l,g12l,g13l,g22l,g23l,g33l, & - 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),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -128,7 +128,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k)) call prim2conM(GRHydro_eos_handle, g11r,g12r,g13r,g22r,g23r,g33r, & - 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),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k), & rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -137,8 +137,8 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) w_lorentzplus(i,j,k)) if(evolve_entropy.ne.0) then - entropyconsminus(i,j,k) = sqrt(avg_detl)*entropyminus(i,j,k)*w_lorentzminus(i,j,k) - entropyconsplus(i,j,k) = sqrt(avg_detr)*entropyplus(i,j,k)*w_lorentzplus(i,j,k) + entropyconsminus(i,j,k) = avg_sdetl*entropyminus(i,j,k)*w_lorentzminus(i,j,k) + entropyconsplus(i,j,k) = avg_sdetr*entropyplus(i,j,k)*w_lorentzplus(i,j,k) end if end do @@ -148,7 +148,7 @@ subroutine primitive2conservativeM(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 g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) @@ -168,8 +168,8 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l,g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r,g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l,g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r,g23r,g33r)) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -187,7 +187,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) endif call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), g11l,g12l,g13l,g22l,g23l,g33l, & - 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),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& @@ -203,7 +203,7 @@ subroutine primitive2conservativeM(CCTK_ARGUMENTS) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), g11r,g12r,g13r,g22r,g23r,g33r, & - 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),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& @@ -235,15 +235,17 @@ end subroutine primitive2conservativeM @@*/ -subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & - dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) +subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, & + x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, 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, & dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz,& @@ -362,18 +364,18 @@ subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy 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 - dBconsx = sqrt(det)*dBvcx - dBconsy = sqrt(det)*dBvcy - dBconsz = sqrt(det)*dBvcz + dBconsx = sdet*dBvcx + dBconsy = sdet*dBvcy + dBconsz = sdet*dBvcz !!$OMP CRITICAL !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w @@ -382,15 +384,16 @@ subroutine prim2conM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy end subroutine prim2conM_hot -subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & - dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) +subroutine prim2conM(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdet, ddens, & + dsx, dsy, dsz, dtau , dBconsx, dBconsy, dBconsz, drho, & + dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) 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, & dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz,& @@ -437,18 +440,18 @@ subroutine prim2conM(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 - dBconsx = sqrt(det)*dBvcx - dBconsy = sqrt(det)*dBvcy - dBconsz = sqrt(det)*dBvcz + dBconsx = sdet*dBvcx + dBconsy = sdet*dBvcy + dBconsz = sdet*dBvcz end subroutine prim2conM @@ -484,7 +487,7 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet CCTK_REAL :: maxtau0 ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -517,25 +520,24 @@ subroutine Primitive2ConservativeCellsM(CCTK_ARGUMENTS) 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 = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ - g22(i,j,k),g23(i,j,k),g33(i,j,k)) + sdet = sdetg(i,j,k) call prim2conM(GRHydro_eos_handle,g11(i,j,k),& g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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), Bvecz(i,j,k), w_lorentz(i,j,k)) if(evolve_entropy.ne.0) then - entropycons(i,j,k) = sqrt(det)*entropy(i,j,k)*w_lorentz(i,j,k) + entropycons(i,j,k) = sdet*entropy(i,j,k)*w_lorentz(i,j,k) end if maxtau0 = max(maxtau0,tau(i,j,k)) @@ -558,19 +560,18 @@ subroutine Primitive2ConservativeCellsM(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 = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ - g22(i,j,k),g23(i,j,k),g33(i,j,k)) + sdet = sdetg(i,j,k) call prim2conM_hot(GRHydro_eos_handle,GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& g11(i,j,k),g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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), & @@ -630,8 +631,8 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -661,8 +662,8 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) ! 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, g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r,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) @@ -680,12 +681,12 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2conpolytypeM(GRHydro_eos_handle, g11l,g12l,g13l,& g22l,g23l,g33l, & - 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),& Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k),& rhominus(i,j,k),& @@ -695,7 +696,7 @@ subroutine Prim2ConservativePolytypeM(CCTK_ARGUMENTS) call prim2conpolytypeM(GRHydro_eos_handle, g11r,g12r,g13r,& g22r,g23r,g33r, & - 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),& Bconsxplus(i,j,k),Bconsyplus(i,j,k),Bconszplus(i,j,k),& rhoplus(i,j,k),& @@ -730,14 +731,14 @@ end subroutine Prim2ConservativePolytypeM @@*/ subroutine prim2conpolytypeM(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & + gzz, sdet, ddens, dsx, dsy, dsz, dtau, dBconsx, dBconsy, dBconsz, & 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, dBconsx, dBconsy, dBconsz, & drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & w_tmp, w, vlowx, vlowy, vlowz @@ -790,18 +791,18 @@ subroutine prim2conpolytypeM(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 - dBconsx = sqrt(det)*dBvcx - dBconsy = sqrt(det)*dBvcy - dBconsz = sqrt(det)*dBvcz + dBconsx = sdet*dBvcx + dBconsy = sdet*dBvcy + dBconsz = sdet*dBvcz end subroutine prim2conpolytypeM @@ -838,7 +839,7 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: sdet ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -865,17 +866,16 @@ subroutine Primitive2ConservativePolyCellsM(CCTK_ARGUMENTS) Bprim => Bvec end if - !$OMP PARALLEL DO PRIVATE(k,j,i,det) + !$OMP PARALLEL DO PRIVATE(k,j,i,sdet) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ - g22(i,j,k),g23(i,j,k),g33(i,j,k)) + sdet = sdetg(i,j,k) call prim2conpolytypeM(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k),& g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(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),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(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), & |