[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 259)
cott at tapir.caltech.edu
cott at tapir.caltech.edu
Tue Aug 2 16:32:01 CDT 2011
User: cott
Date: 2011/08/02 04:32 PM
Modified:
/trunk/src/
GRHydro_FluxSplit.F90, GRHydro_HLLE.F90, GRHydro_HLLEM.F90, GRHydro_Marquina.F90, GRHydro_ParamCheck.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructM.F90, GRHydro_ReconstructPoly.F90, GRHydro_ReconstructPolyM.F90, GRHydro_RoeSolver.F90, GRHydro_Source.F90, GRHydro_SourceM.F90, GRHydro_Tmunu.F90, GRHydro_TmunuM.F90
Log:
* Optimize: remove support for shift_state = 0 (except for shock tubes and
Cowling calculations of spherically symmetric objects, there is no reason
not to have storage for the shift).
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_FluxSplit.F90
Delta lines: +9 -31
===================================================================
--- trunk/src/GRHydro_FluxSplit.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_FluxSplit.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -73,11 +73,7 @@
gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
gyz(i,j,k),gzz(i,j,k))
- if (shift_state .eq. 0) then
- beta = 0.d0
- else
- beta = betax(i,j,k)
- end if
+ beta = betax(i,j,k)
call eigenvalues(GRHydro_eos_handle, &
rho (i,j,k), &
@@ -119,11 +115,7 @@
gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
gyz(i,j,k),gzz(i,j,k))
- if (shift_state .eq. 0) then
- beta = 0.d0
- else
- beta = betay(i,j,k)
- end if
+ beta = betay(i,j,k)
call eigenvalues(GRHydro_eos_handle, &
rho (i,j,k), &
@@ -165,11 +157,7 @@
gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),&
gyz(i,j,k),gzz(i,j,k))
- if (shift_state .eq. 0) then
- beta = 0.d0
- else
- beta = betaz(i,j,k)
- end if
+ beta = betaz(i,j,k)
call eigenvalues(GRHydro_eos_handle, &
rho (i,j,k), &
@@ -310,11 +298,8 @@
do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 2
do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 2
- if (shift_state .eq. 0) then
- dummy = 0.d0
- else
- dummy = betax(:,j,k)
- end if
+ dummy = betax(:,j,k)
+
do i = 1, cctk_lsh(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))
@@ -347,11 +332,8 @@
do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 2
do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 2
- if (shift_state .eq. 0) then
- dummy = 0.d0
- else
- dummy = betay(i,:,k)
- end if
+ dummy = betay(i,:,k)
+
do j = 1, cctk_lsh(2)
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))
@@ -384,12 +366,8 @@
do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 2
do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 2
- if (shift_state .eq. 0) then
- dummy = 0.d0
- else
- dummy = betaz(i,j,:)
- end if
-
+ dummy = betaz(i,j,:)
+
do k = 1, cctk_lsh(3)
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))
File [modified]: GRHydro_HLLE.F90
Delta lines: +21 -29
===================================================================
--- trunk/src/GRHydro_HLLE.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_HLLE.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -103,22 +103,18 @@
!!$ Note also need the average of the lapse at the
!!$ left and right points.
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
+ if (flux_direction == 1) then
+ avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+ betax(i,j,k))
+ else if (flux_direction == 2) then
+ avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+ betay(i,j,k))
+ else if (flux_direction == 3) then
+ avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+ betaz(i,j,k))
else
- avg_beta = 0.d0
- end if
+ call CCTK_WARN(0, "Flux direction not x,y,z")
+ end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
@@ -470,21 +466,17 @@
!!$ Note also need the average of the lapse at the
!!$ left and right points.
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
+ if (flux_direction == 1) then
+ avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+ betax(i,j,k))
+ else if (flux_direction == 2) then
+ avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+ betay(i,j,k))
+ else if (flux_direction == 3) then
+ avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+ betaz(i,j,k))
else
- avg_beta = 0.d0
+ call CCTK_WARN(0, "Flux direction not x,y,z")
end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
File [modified]: GRHydro_HLLEM.F90
Delta lines: +23 -34
===================================================================
--- trunk/src/GRHydro_HLLEM.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_HLLEM.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -144,27 +144,20 @@
!!$
!!$ In MHD, we need all three shift components regardless of the flux direction
- if (shift_state .ne. 0) then
- avg_betax = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- avg_betay = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- avg_betaz = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- if (flux_direction == 1) then
- avg_beta=avg_betax
- else if (flux_direction == 2) then
- avg_beta=avg_betay
- else if (flux_direction == 3) then
- avg_beta=avg_betaz
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
+ avg_betax = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+ betax(i,j,k))
+ avg_betay = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+ betay(i,j,k))
+ avg_betaz = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+ betaz(i,j,k))
+ if (flux_direction == 1) then
+ avg_beta=avg_betax
+ else if (flux_direction == 2) then
+ avg_beta=avg_betay
+ else if (flux_direction == 3) then
+ avg_beta=avg_betaz
else
- avg_beta = 0.d0
- avg_betax = 0.d0
- avg_betay = 0.d0
- avg_betaz = 0.d0
+ call CCTK_WARN(0, "Flux direction not x,y,z")
end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
@@ -580,21 +573,17 @@
!!$ Note also need the average of the lapse at the
!!$ left and right points.
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
+ if (flux_direction == 1) then
+ avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+ betax(i,j,k))
+ else if (flux_direction == 2) then
+ avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+ betay(i,j,k))
+ else if (flux_direction == 3) then
+ avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+ betaz(i,j,k))
else
- avg_beta = 0.d0
+ call CCTK_WARN(0, "Flux direction not x,y,z")
end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
File [modified]: GRHydro_Marquina.F90
Delta lines: +8 -12
===================================================================
--- trunk/src/GRHydro_Marquina.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_Marquina.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -111,19 +111,15 @@
!!$ Set metric terms at interface
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
- betay(i,j,k))
- else
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- end if
+ if (flux_direction == 1) then
+ avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+ betax(i,j,k))
+ else if (flux_direction == 2) then
+ avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+ betay(i,j,k))
else
- avg_beta = 0.d0
+ avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+ betaz(i,j,k))
end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
File [modified]: GRHydro_ParamCheck.F90
Delta lines: +4 -0
===================================================================
--- trunk/src/GRHydro_ParamCheck.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_ParamCheck.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -92,6 +92,10 @@
evolve_MHD = 0
endif
+ if(shift_state.eq.0) then
+ call CCTK_PARAMWARN("shift_state = 0 (no shift storage) no longer supported!d");
+ endif
+
if (CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then
evolve_Y_e = 1
else
File [modified]: GRHydro_Reconstruct.F90
Delta lines: +0 -5
===================================================================
--- trunk/src/GRHydro_Reconstruct.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_Reconstruct.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -41,11 +41,6 @@
CCTK_INT :: i,j,k
CCTK_REAL :: local_min_tracer
- ! set things up
- if (shift_state .eq. 0) then
- call CCTK_WARN(0,"This code no longer supports shift_state = 0");
- endif
-
if (CCTK_EQUALS(recon_method,"tvd")) then
! this handles MHD and non-MHD
call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF)
File [modified]: GRHydro_ReconstructM.F90
Delta lines: +3 -9
===================================================================
--- trunk/src/GRHydro_ReconstructM.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_ReconstructM.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -107,15 +107,9 @@
psi4 = 1.d0
- if (shift_state .ne. 0) then
- lbetax = betax
- lbetay = betay
- lbetaz = betaz
- else
- lbetax = 0.d0
- lbetay = 0.d0
- lbetaz = 0.d0
- end if
+ lbetax = betax
+ lbetay = betay
+ lbetaz = betaz
!!$ Initialize variables that store reconstructed quantities
File [modified]: GRHydro_ReconstructPoly.F90
Delta lines: +6 -13
===================================================================
--- trunk/src/GRHydro_ReconstructPoly.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_ReconstructPoly.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -118,20 +118,13 @@
!$OMP WORKSHARE
psi4 = 1.0d0
!$OMP END WORKSHARE NOWAIT
- if (shift_state .ne. 0) then
- !$OMP WORKSHARE
- lbetax = betax
- lbetay = betay
- lbetaz = betaz
- !$OMP END WORKSHARE NOWAIT
- else
- !$OMP WORKSHARE
- lbetax = 0.d0
- lbetay = 0.d0
- lbetaz = 0.d0
- !$OMP END WORKSHARE NOWAIT
- end if
+ !$OMP WORKSHARE
+ lbetax = betax
+ lbetay = betay
+ lbetaz = betaz
+ !$OMP END WORKSHARE NOWAIT
+
!!$ Initialize variables that store reconstructed quantities
!$OMP WORKSHARE
File [modified]: GRHydro_ReconstructPolyM.F90
Delta lines: +3 -9
===================================================================
--- trunk/src/GRHydro_ReconstructPolyM.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_ReconstructPolyM.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -107,15 +107,9 @@
psi4 = 1.d0
- if (shift_state .ne. 0) then
- lbetax = betax
- lbetay = betay
- lbetaz = betaz
- else
- lbetax = 0.d0
- lbetay = 0.d0
- lbetaz = 0.d0
- end if
+ lbetax = betax
+ lbetay = betay
+ lbetaz = betaz
!!$ Initialize variables that store reconstructed quantities
File [modified]: GRHydro_RoeSolver.F90
Delta lines: +9 -13
===================================================================
--- trunk/src/GRHydro_RoeSolver.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_RoeSolver.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -107,21 +107,17 @@
!!$ Set metric terms at interface
- if (shift_state .ne. 0) then
- if (flux_direction == 1) then
- avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
- betax(i,j,k))
- else if (flux_direction == 2) then
- avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
+ if (flux_direction == 1) then
+ avg_beta = 0.5d0 * (betax(i+xoffset,j+yoffset,k+zoffset) + &
+ betax(i,j,k))
+ else if (flux_direction == 2) then
+ avg_beta = 0.5d0 * (betay(i+xoffset,j+yoffset,k+zoffset) + &
betay(i,j,k))
- else if (flux_direction == 3) then
- avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
- betaz(i,j,k))
- else
- call CCTK_WARN(0, "Flux direction not x,y,z")
- end if
+ else if (flux_direction == 3) then
+ avg_beta = 0.5d0 * (betaz(i+xoffset,j+yoffset,k+zoffset) + &
+ betaz(i,j,k))
else
- avg_beta = 0.d0
+ call CCTK_WARN(0, "Flux direction not x,y,z")
end if
avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset))
File [modified]: GRHydro_Source.F90
Delta lines: +30 -52
===================================================================
--- trunk/src/GRHydro_Source.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_Source.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -175,62 +175,40 @@
rhoenthalpyW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k))*&
w_lorentz(i,j,k)**2
- if (shift_state .ne. 0) then
+ shiftx = betax(i,j,k)
+ shifty = betay(i,j,k)
+ shiftz = betaz(i,j,k)
- shiftx = betax(i,j,k)
- shifty = betay(i,j,k)
- shiftz = betaz(i,j,k)
+ if (local_spatial_order .eq. 2) then
+
+ dx_betax = DIFF_X_2(betax)
+ dx_betay = DIFF_X_2(betay)
+ dx_betaz = DIFF_X_2(betaz)
+
+ dy_betax = DIFF_Y_2(betax)
+ dy_betay = DIFF_Y_2(betay)
+ dy_betaz = DIFF_Y_2(betaz)
+
+ dz_betax = DIFF_Z_2(betax)
+ dz_betay = DIFF_Z_2(betay)
+ dz_betaz = DIFF_Z_2(betaz)
+
+ else
- if (local_spatial_order .eq. 2) then
-
- dx_betax = DIFF_X_2(betax)
- dx_betay = DIFF_X_2(betay)
- dx_betaz = DIFF_X_2(betaz)
+ dx_betax = DIFF_X_4(betax)
+ dx_betay = DIFF_X_4(betay)
+ dx_betaz = DIFF_X_4(betaz)
+
+ dy_betax = DIFF_Y_4(betax)
+ dy_betay = DIFF_Y_4(betay)
+ dy_betaz = DIFF_Y_4(betaz)
- dy_betax = DIFF_Y_2(betax)
- dy_betay = DIFF_Y_2(betay)
- dy_betaz = DIFF_Y_2(betaz)
-
- dz_betax = DIFF_Z_2(betax)
- dz_betay = DIFF_Z_2(betay)
- dz_betaz = DIFF_Z_2(betaz)
-
- else
-
- dx_betax = DIFF_X_4(betax)
- dx_betay = DIFF_X_4(betay)
- dx_betaz = DIFF_X_4(betaz)
-
- dy_betax = DIFF_Y_4(betax)
- dy_betay = DIFF_Y_4(betay)
- dy_betaz = DIFF_Y_4(betaz)
-
- dz_betax = DIFF_Z_4(betax)
- dz_betay = DIFF_Z_4(betay)
- dz_betaz = DIFF_Z_4(betaz)
-
- end if
+ dz_betax = DIFF_Z_4(betax)
+ dz_betay = DIFF_Z_4(betay)
+ dz_betaz = DIFF_Z_4(betaz)
+
+ end if
- else
-
- shiftx = 0.0d0
- shifty = 0.0d0
- shiftz = 0.0d0
-
- dx_betax = 0.0d0
- dx_betay = 0.0d0
- dx_betaz = 0.0d0
-
- dy_betax = 0.0d0
- dy_betay = 0.0d0
- dy_betaz = 0.0d0
-
- dz_betax = 0.0d0
- dz_betay = 0.0d0
- dz_betaz = 0.0d0
-
- endif
-
invalp = 1.0d0 / alp(i,j,k)
invalp2 = invalp**2
velxshift = velx(i,j,k) - shiftx*invalp
File [modified]: GRHydro_SourceM.F90
Delta lines: +29 -51
===================================================================
--- trunk/src/GRHydro_SourceM.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_SourceM.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -199,62 +199,40 @@
localgxy, localgxz, localgyy, localgyz, localgzz)
- if (shift_state .ne. 0) then
-
- shiftx = betax(i,j,k)
- shifty = betay(i,j,k)
- shiftz = betaz(i,j,k)
-
- if (local_spatial_order .eq. 2) then
-
- dx_betax = DIFF_X_2(betax)
- dx_betay = DIFF_X_2(betay)
- dx_betaz = DIFF_X_2(betaz)
+ shiftx = betax(i,j,k)
+ shifty = betay(i,j,k)
+ shiftz = betaz(i,j,k)
+
+ if (local_spatial_order .eq. 2) then
+
+ dx_betax = DIFF_X_2(betax)
+ dx_betay = DIFF_X_2(betay)
+ dx_betaz = DIFF_X_2(betaz)
+
+ dy_betax = DIFF_Y_2(betax)
+ dy_betay = DIFF_Y_2(betay)
+ dy_betaz = DIFF_Y_2(betaz)
- dy_betax = DIFF_Y_2(betax)
- dy_betay = DIFF_Y_2(betay)
- dy_betaz = DIFF_Y_2(betaz)
-
- dz_betax = DIFF_Z_2(betax)
- dz_betay = DIFF_Z_2(betay)
- dz_betaz = DIFF_Z_2(betaz)
+ dz_betax = DIFF_Z_2(betax)
+ dz_betay = DIFF_Z_2(betay)
+ dz_betaz = DIFF_Z_2(betaz)
- else
-
- dx_betax = DIFF_X_4(betax)
- dx_betay = DIFF_X_4(betay)
- dx_betaz = DIFF_X_4(betaz)
-
- dy_betax = DIFF_Y_4(betax)
- dy_betay = DIFF_Y_4(betay)
- dy_betaz = DIFF_Y_4(betaz)
-
- dz_betax = DIFF_Z_4(betax)
- dz_betay = DIFF_Z_4(betay)
- dz_betaz = DIFF_Z_4(betaz)
-
- end if
-
else
- shiftx = 0.0d0
- shifty = 0.0d0
- shiftz = 0.0d0
-
- dx_betax = 0.0d0
- dx_betay = 0.0d0
- dx_betaz = 0.0d0
+ dx_betax = DIFF_X_4(betax)
+ dx_betay = DIFF_X_4(betay)
+ dx_betaz = DIFF_X_4(betaz)
+
+ dy_betax = DIFF_Y_4(betax)
+ dy_betay = DIFF_Y_4(betay)
+ dy_betaz = DIFF_Y_4(betaz)
+
+ dz_betax = DIFF_Z_4(betax)
+ dz_betay = DIFF_Z_4(betay)
+ dz_betaz = DIFF_Z_4(betaz)
+
+ end if
- dy_betax = 0.0d0
- dy_betay = 0.0d0
- dy_betaz = 0.0d0
-
- dz_betax = 0.0d0
- dz_betay = 0.0d0
- dz_betaz = 0.0d0
-
- endif
-
invalp = 1.0d0 / alp(i,j,k)
invalp2 = invalp**2
velxshift = velx(i,j,k) - shiftx*invalp
File [modified]: GRHydro_Tmunu.F90
Delta lines: +5 -17
===================================================================
--- trunk/src/GRHydro_Tmunu.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_Tmunu.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -92,24 +92,12 @@
!!$ Calculate lower components and square of shift vector.
- if (shift_state .ne. 0) then
-
- betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
- betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
- betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
-
- beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow
+ betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
+ betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
+ betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
+
+ beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow
- else
-
- betaxlow = 0.0D0
- betaylow = 0.0D0
- betazlow = 0.0D0
-
- beta2 = 0.0D0
-
- end if
-
!!$ Calculate the specific relativistic enthalpy times rho times the
!!$ square of the lorentz factor.
File [modified]: GRHydro_TmunuM.F90
Delta lines: +12 -27
===================================================================
--- trunk/src/GRHydro_TmunuM.F90 2011-08-02 20:19:41 UTC (rev 258)
+++ trunk/src/GRHydro_TmunuM.F90 2011-08-02 21:32:00 UTC (rev 259)
@@ -100,36 +100,21 @@
!!$ Calculate lower components and square of shift vector.
- if (shift_state .ne. 0) then
-
- betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
- betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
- betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
- beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow
-
- bdotbeta = betaxlow*Bvecx(i,j,k)+betaylow*Bvecy(i,j,k)+betazlow*Bvecz(i,j,k)
- vdotbeta = betaxlow*velx(i,j,k)+betaylow*vely(i,j,k)+betazlow*velz(i,j,k)
+
+ betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k)
+ betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k)
+ betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k)
+ beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow
+
+ bdotbeta = betaxlow*Bvecx(i,j,k)+betaylow*Bvecy(i,j,k)+betazlow*Bvecz(i,j,k)
+ vdotbeta = betaxlow*velx(i,j,k)+betaylow*vely(i,j,k)+betazlow*velz(i,j,k)
!!$ u0 low is missing the w_lorentz factor (see below)!!
- utlow = -1.d0*alp(i,j,k) + vdotbeta
+ utlow = -1.d0*alp(i,j,k) + vdotbeta
+
+ btlow = -1.0d0*w_lorentz(i,j,k)*alp(i,j,k)*bdotv + &
+ bdotbeta/w_lorentz(i,j,k) + w_lorentz(i,j,k)*bdotv*vdotbeta
- btlow = -1.0d0*w_lorentz(i,j,k)*alp(i,j,k)*bdotv + &
- bdotbeta/w_lorentz(i,j,k) + w_lorentz(i,j,k)*bdotv*vdotbeta
-
-
- else
-
- betaxlow = 0.0D0
- betaylow = 0.0D0
- betazlow = 0.0D0
- beta2 = 0.0D0
-
-!!$ u0 low is missing the w_lorentz factor (see below)!!
- utlow = -1.0*alp(i,j,k)
- btlow = utlow*w_lorentz(i,j,k)*bdotv
-
- end if
-
!!$ Calculate the specific relativistic enthalpy times rho + the mag. field contribution times the
!!$ square of the lorentz factor.
More information about the Commits
mailing list