diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-01-14 14:23:33 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-01-14 14:23:33 +0000 |
| commit | 48b053d41b2aa9ce720d443d288eef535fb9651d (patch) | |
| tree | bb5274de4dfb94de758bafc0ee77e6b5f499fb06 /src/GRHydro_Reconstruct.F90 | |
| parent | d6e45e79d92f3d2adf3eb4ebcc7641db07119f30 (diff) | |
GRHydro: Implement entropy evolution (separated from temp evolution)
From: Bruno Coutinho Mundim <bcmsma@astro.rit.edu>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@455 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat (limited to 'src/GRHydro_Reconstruct.F90')
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 62 |
1 files changed, 31 insertions, 31 deletions
diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index 7168523..8ac033a 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -132,43 +132,43 @@ subroutine Reconstruction(CCTK_ARGUMENTS) epsminus(i,j,k) = eps(i,j,k) if (reconstruct_Wv.ne.0) then ! divide out the Loretnz factor for both the - ! plus and minus quantities this should by construction ensure - ! that any Lorentz factor calculated from them later on is - ! physical (ie. > 1.d0) - velxplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1) + ! plus and minus quantities this should by construction ensure + ! that any Lorentz factor calculated from them later on is + ! physical (ie. > 1.d0) + velxplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1) velyplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2) velzplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3) agxx = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) - agxy = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) - agxz = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) - agyy = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) - agyz = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) - agzz = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - w = sqrt( 1.d0 + agxx*velxplus(i,j,k)*velxplus(i,j,k) + & - agyy*velyplus(i,j,k)*velyplus(i,j,k) & - + agzz*velzplus(i,j,k)*velzplus(i,j,k) + & - 2.d0*agxy*velxplus(i,j,k)*velyplus(i,j,k) & - + 2.d0*agxz*velxplus(i,j,k)*velzplus(i,j,k) + & - 2.d0*agyz*velyplus(i,j,k)*velzplus(i,j,k) ) - velxplus(i,j,k) = velxplus(i,j,k)/w - velyplus(i,j,k) = velyplus(i,j,k)/w - velzplus(i,j,k) = velzplus(i,j,k)/w - + agxy = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + agxz = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + agyy = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + agyz = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + agzz = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) + w = sqrt( 1.d0 + agxx*velxplus(i,j,k)*velxplus(i,j,k) & + + agyy*velyplus(i,j,k)*velyplus(i,j,k) & + + agzz*velzplus(i,j,k)*velzplus(i,j,k) & + + 2.d0*agxy*velxplus(i,j,k)*velyplus(i,j,k) & + + 2.d0*agxz*velxplus(i,j,k)*velzplus(i,j,k) & + + 2.d0*agyz*velyplus(i,j,k)*velzplus(i,j,k) ) + velxplus(i,j,k) = velxplus(i,j,k)/w + velyplus(i,j,k) = velyplus(i,j,k)/w + velzplus(i,j,k) = velzplus(i,j,k)/w + velxminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1) velyminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2) velzminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3) - agxx = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) - agxy = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) - agxz = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) - agyy = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) - agyz = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) - agzz = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) - w = sqrt( 1.d0 + agxx*velxminus(i,j,k)*velxminus(i,j,k) + & - agyy*velyminus(i,j,k)*velyminus(i,j,k) & - + agzz*velzminus(i,j,k)*velzminus(i,j,k) + & - 2.d0*agxy*velxminus(i,j,k)*velyminus(i,j,k) & - + 2.d0*agxz*velxminus(i,j,k)*velzminus(i,j,k) + & - 2.d0*agyz*velyminus(i,j,k)*velzminus(i,j,k) ) + agxx = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + agxy = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + agxz = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + agyy = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + agyz = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + agzz = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + w = sqrt( 1.d0 + agxx*velxminus(i,j,k)*velxminus(i,j,k) & + + agyy*velyminus(i,j,k)*velyminus(i,j,k) & + + agzz*velzminus(i,j,k)*velzminus(i,j,k) & + + 2.d0*agxy*velxminus(i,j,k)*velyminus(i,j,k) & + + 2.d0*agxz*velxminus(i,j,k)*velzminus(i,j,k) & + + 2.d0*agyz*velyminus(i,j,k)*velzminus(i,j,k) ) velxminus(i,j,k) = velxminus(i,j,k)/w velyminus(i,j,k) = velyminus(i,j,k)/w velzminus(i,j,k) = velzminus(i,j,k)/w |