diff options
Diffstat (limited to 'src/GRHydro_Source.F90')
| -rw-r--r-- | src/GRHydro_Source.F90 | 82 |
1 files changed, 30 insertions, 52 deletions
diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90 index 1a981db..f8fd9eb 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -175,62 +175,40 @@ subroutine SourceTerms(CCTK_ARGUMENTS) 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 - - 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 - + 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 - 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 + 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 - endif - invalp = 1.0d0 / alp(i,j,k) invalp2 = invalp**2 velxshift = velx(i,j,k) - shiftx*invalp |