c Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations c Copyright (C) 2001 Mark Miller /*@@ @file hydro_source_update.F @date December 1999 @author Mark Miller @desc Calculate the source terms for the hydro update. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "CactusEinstein/Einstein/src/Einstein.h" /*@@ @routine hydro_source_update @date December 1999 @author Mark Miller @desc Calculate the source terms for the hydro update. @enddesc @calls @calledby @history @endhistory @@*/ subroutine hydro_source_update(CCTK_FARGUMENTS,evolve_dt) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_REAL evolve_dt c LOCAL VARS integer i,j,k,nx,ny,nz CCTK_REAL one,two,half CCTK_REAL t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz CCTK_REAL chris_xx0,chris_xy0,chris_xz0,chris_yy0,chris_yz0,chris_zz0 CCTK_REAL sqrtdet,det,uxx,uxy,uxz,uyy,uyz,uzz,rhoenthalpy CCTK_REAL shiftx,shifty,shiftz,velxshift,velyshift,velzshift CCTK_REAL vlowx,vlowy,vlowz CCTK_REAL dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz CCTK_REAL dz_betax,dz_betay,dz_betaz CCTK_REAL dx_alp,dy_alp,dz_alp CCTK_REAL tau_source,sx_source,sy_source,sz_source CCTK_REAL dx_gxx,dx_gxy,dx_gxz,dx_gyy,dx_gyz,dx_gzz CCTK_REAL dy_gxx,dy_gxy,dy_gxz,dy_gyy,dy_gyz,dy_gzz CCTK_REAL dz_gxx,dz_gxy,dz_gxz,dz_gyy,dz_gyz,dz_gzz CCTK_REAL dx,dy,dz one = 1.0d0 two = 2.0d0 half = 0.5d0 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 c calculate some aux variables det= -(gxz(i,j,k)**2*gyy(i,j,k)) + . two*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) - . gxx(i,j,k)*gyz(i,j,k)**2 . - gxy(i,j,k)**2*gzz(i,j,k) + . gxx(i,j,k)*gyy(i,j,k)*gzz(i,j,k) sqrtdet = sqrt(det) uxx=(-gyz(i,j,k)**2 + gyy(i,j,k)*gzz(i,j,k))/det uxy=(gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k))/det uyy=(-gxz(i,j,k)**2 + gxx(i,j,k)*gzz(i,j,k))/det uxz=(-gxz(i,j,k)*gyy(i,j,k) + gxy(i,j,k)*gyz(i,j,k))/det uyz=(gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k))/det uzz=(-gxy(i,j,k)**2 + gxx(i,j,k)*gyy(i,j,k))/det rhoenthalpy = rho(i,j,k) * (one + eps(i,j,k)) + press(i,j,k) if(shift_state .eq. SHIFT_INACTIVE) then shiftx = 0.0d0 shifty = 0.0d0 shiftz = 0.0d0 c shift derivatives 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 else shiftx = betax(i,j,k) shifty = betay(i,j,k) shiftz = betaz(i,j,k) c shift derivatives dx_betax = (betax(i+1,j,k) - betax(i-1,j,k))/(2.0D0*dx) dx_betay = (betay(i+1,j,k) - betay(i-1,j,k))/(2.0D0*dx) dx_betaz = (betaz(i+1,j,k) - betaz(i-1,j,k))/(2.0D0*dx) dy_betax = (betax(i,j+1,k) - betax(i,j-1,k))/(2.0D0*dy) dy_betay = (betay(i,j+1,k) - betay(i,j-1,k))/(2.0D0*dy) dy_betaz = (betaz(i,j+1,k) - betaz(i,j-1,k))/(2.0D0*dy) dz_betax = (betax(i,j,k+1) - betax(i,j,k-1))/(2.0D0*dz) dz_betay = (betay(i,j,k+1) - betay(i,j,k-1))/(2.0D0*dz) dz_betaz = (betaz(i,j,k+1) - betaz(i,j,k-1))/(2.0D0*dz) endif velxshift = velx(i,j,k) - shiftx/alp(i,j,k) velyshift = vely(i,j,k) - shifty/alp(i,j,k) velzshift = velz(i,j,k) - shiftz/alp(i,j,k) vlowx = velx(i,j,k)*gxx(i,j,k) + vely(i,j,k)*gxy(i,j,k) + . velz(i,j,k)*gxz(i,j,k) vlowy = velx(i,j,k)*gxy(i,j,k) + vely(i,j,k)*gyy(i,j,k) + . velz(i,j,k)*gyz(i,j,k) vlowz = velx(i,j,k)*gxz(i,j,k) + vely(i,j,k)*gyz(i,j,k) + . velz(i,j,k)*gzz(i,j,k) c calculate the stress energy tensor.... c (these indices are up) t00 = ( rhoenthalpy * w_lorentz(i,j,k)**2 - press(i,j,k))/ . (alp(i,j,k)**2) t0x = rhoenthalpy*w_lorentz(i,j,k)**2 * velxshift/alp(i,j,k)+ . press(i,j,k) * shiftx / (alp(i,j,k)**2) t0y = rhoenthalpy*w_lorentz(i,j,k)**2 * velyshift/alp(i,j,k)+ . press(i,j,k) * shifty / (alp(i,j,k)**2) t0z = rhoenthalpy*w_lorentz(i,j,k)**2 * velzshift/alp(i,j,k)+ . press(i,j,k) * shiftz / (alp(i,j,k)**2) txx = rhoenthalpy*w_lorentz(i,j,k)**2*velxshift*velxshift + . press(i,j,k) * (uxx - shiftx*shiftx/(alp(i,j,k)**2)) txy = rhoenthalpy*w_lorentz(i,j,k)**2*velxshift*velyshift + . press(i,j,k) * (uxy - shiftx*shifty/(alp(i,j,k)**2)) txz = rhoenthalpy*w_lorentz(i,j,k)**2*velxshift*velzshift + . press(i,j,k) * (uxz - shiftx*shiftz/(alp(i,j,k)**2)) tyy = rhoenthalpy*w_lorentz(i,j,k)**2*velyshift*velyshift + . press(i,j,k) * (uyy - shifty*shifty/(alp(i,j,k)**2)) tyz = rhoenthalpy*w_lorentz(i,j,k)**2*velyshift*velzshift + . press(i,j,k) * (uyz - shifty*shiftz/(alp(i,j,k)**2)) tzz = rhoenthalpy*w_lorentz(i,j,k)**2*velzshift*velzshift + . press(i,j,k) * (uzz - shiftz*shiftz/(alp(i,j,k)**2)) c lapse derivatives dx_alp = (alp(i+1,j,k) - alp(i-1,j,k))/(2.0D0*dx) dy_alp = (alp(i,j+1,k) - alp(i,j-1,k))/(2.0D0*dy) dz_alp = (alp(i,j,k+1) - alp(i,j,k-1))/(2.0D0*dz) c metric derivatives dx_gxx = (gxx(i+1,j,k) - gxx(i-1,j,k))/(2.0D0*dx) dx_gxy = (gxy(i+1,j,k) - gxy(i-1,j,k))/(2.0D0*dx) dx_gxz = (gxz(i+1,j,k) - gxz(i-1,j,k))/(2.0D0*dx) dx_gyy = (gyy(i+1,j,k) - gyy(i-1,j,k))/(2.0D0*dx) dx_gyz = (gyz(i+1,j,k) - gyz(i-1,j,k))/(2.0D0*dx) dx_gzz = (gzz(i+1,j,k) - gzz(i-1,j,k))/(2.0D0*dx) dy_gxx = (gxx(i,j+1,k) - gxx(i,j-1,k))/(2.0D0*dy) dy_gxy = (gxy(i,j+1,k) - gxy(i,j-1,k))/(2.0D0*dy) dy_gxz = (gxz(i,j+1,k) - gxz(i,j-1,k))/(2.0D0*dy) dy_gyy = (gyy(i,j+1,k) - gyy(i,j-1,k))/(2.0D0*dy) dy_gyz = (gyz(i,j+1,k) - gyz(i,j-1,k))/(2.0D0*dy) dy_gzz = (gzz(i,j+1,k) - gzz(i,j-1,k))/(2.0D0*dy) dz_gxx = (gxx(i,j,k+1) - gxx(i,j,k-1))/(2.0D0*dz) dz_gxy = (gxy(i,j,k+1) - gxy(i,j,k-1))/(2.0D0*dz) dz_gxz = (gxz(i,j,k+1) - gxz(i,j,k-1))/(2.0D0*dz) dz_gyy = (gyy(i,j,k+1) - gyy(i,j,k-1))/(2.0D0*dz) dz_gyz = (gyz(i,j,k+1) - gyz(i,j,k-1))/(2.0D0*dz) dz_gzz = (gzz(i,j,k+1) - gzz(i,j,k-1))/(2.0D0*dz) c space-space-time Christoffels chris_xx0 = - kxx(i,j,k)/alp(i,j,k) chris_xy0 = - kxy(i,j,k)/alp(i,j,k) chris_xz0 = - kxz(i,j,k)/alp(i,j,k) chris_yy0 = - kyy(i,j,k)/alp(i,j,k) chris_yz0 = - kyz(i,j,k)/alp(i,j,k) chris_zz0 = - kzz(i,j,k)/alp(i,j,k) tau_source = t00 * . (shiftx*shiftx*kxx(i,j,k) + shifty*shifty*kyy(i,j,k) + . shiftz*shiftz*kzz(i,j,k) + two*shiftx*shifty*kxy(i,j,k) + . two*shiftx*shiftz*kxz(i,j,k)+two*shifty*shiftz*kyz(i,j,k)- . (shiftx*dx_alp + shifty*dy_alp + shiftz*dz_alp) ) + . t0x * (- dx_alp + two * (shiftx*kxx(i,j,k) + . shifty*kxy(i,j,k) + shiftz*kxz(i,j,k))) + . t0y * (- dy_alp + two * (shiftx*kxy(i,j,k) + . shifty*kyy(i,j,k) + shiftz*kyz(i,j,k))) + . t0z * (- dz_alp + two * (shiftx*kxz(i,j,k) + . shifty*kyz(i,j,k) + shiftz*kzz(i,j,k))) - . alp(i,j,k) * (txx * chris_xx0 + tyy * chris_yy0 + . tzz * chris_zz0 + two * txy * chris_xy0 + . two * txz * chris_xz0 + two * tyz * chris_yz0) sx_source = t00 * . (half*shiftx*shiftx*dx_gxx + half*shifty*shifty*dx_gyy+ . half*shiftz*shiftz*dx_gzz + shiftx*shifty*dx_gxy + . shiftx*shiftz*dx_gxz + shifty*shiftz*dx_gyz - . alp(i,j,k) * dx_alp) + . t0x * (shiftx*dx_gxx + shifty*dx_gxy + shiftz*dx_gxz) + . t0y * (shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) + . t0z * (shiftx*dx_gxz + shifty*dx_gyz + shiftz*dx_gzz) + . (half*txx*dx_gxx + half*tyy*dx_gyy + half*tzz*dx_gzz + . txy*dx_gxy + txz*dx_gxz + tyz*dx_gyz) + . rhoenthalpy * w_lorentz(i,j,k)**2 * . (vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz)/ . alp(i,j,k) sy_source = t00 * . (half*shiftx*shiftx*dy_gxx + half*shifty*shifty*dy_gyy+ . half*shiftz*shiftz*dy_gzz + shiftx*shifty*dy_gxy + . shiftx*shiftz*dy_gxz + shifty*shiftz*dy_gyz - . alp(i,j,k) * dy_alp) + . t0x * (shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) + . t0y * (shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) + . t0z * (shiftx*dy_gxz + shifty*dy_gyz + shiftz*dy_gzz) + . (half*txx*dy_gxx + half*tyy*dy_gyy + half*tzz*dy_gzz + . txy*dy_gxy + txz*dy_gxz + tyz*dy_gyz) + . rhoenthalpy * w_lorentz(i,j,k)**2 * . (vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz)/ . alp(i,j,k) sz_source = t00 * . (half*shiftx*shiftx*dz_gxx + half*shifty*shifty*dz_gyy+ . half*shiftz*shiftz*dz_gzz + shiftx*shifty*dz_gxy + . shiftx*shiftz*dz_gxz + shifty*shiftz*dz_gyz - . alp(i,j,k) * dz_alp) + . t0x * (shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) + . t0y * (shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) + . t0z * (shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) + . (half*txx*dz_gxx + half*tyy*dz_gyy + half*tzz*dz_gzz + . txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz) + . rhoenthalpy * w_lorentz(i,j,k)**2 * . (vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz)/ . alp(i,j,k) c no source for dens c sx source update del_sx(i,j,k) = del_sx(i,j,k) + evolve_dt * ( . alp(i,j,k) * sqrtdet * sx_source) c sy source update del_sy(i,j,k) = del_sy(i,j,k) + evolve_dt * ( . alp(i,j,k) * sqrtdet * sy_source) c sz source update del_sz(i,j,k) = del_sz(i,j,k) + evolve_dt * ( . alp(i,j,k) * sqrtdet * sz_source) c tau source update del_tau(i,j,k) = del_tau(i,j,k) + evolve_dt * ( . alp(i,j,k) * sqrtdet * tau_source) enddo enddo enddo return end