c Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations c Copyright (C) 2001 Mark Miller /*@@ @file roe_flux_x_update.F @date December 1999 @author Mark Miller @desc Performs the X flux update using Roe method. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" #include "CactusEinstein/Einstein/src/Einstein.h" /*@@ @routine roe_flux_x_update @date December 1999 @author Mark Miller @desc can do reconstruction based on primative or evolved variables. For reconstruction based on primates, #define RECON_PRIM For reconstruction based on evolved vars, #define RECON_EVOL @enddesc @calls @calledby @history @endhistory @@*/ #define RECON_PRIM subroutine roe_flux_x_update(CCTK_FARGUMENTS,evolve_dt,order) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS CCTK_REAL evolve_dt integer order c LOCAL VARS CCTK_REAL num_dens_flux(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL num_tau_flux(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL num_sx_flux(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL num_sy_flux(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) CCTK_REAL num_sz_flux(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) integer i,j,k,tvd,fmi,nx,ny,nz,CCTK_Equals logical eno CCTK_REAL dpdrho,dpdeps CCTK_REAL betaxh,betayh,betazh,bxxh,bxyh,bxzh,byxh,byyh,byzh CCTK_REAL bzxh,bzyh,bzzh CCTK_REAL slopel,sloper CCTK_REAL dens_r,tau_r,sx_r,sy_r,sz_r,rho_r,eps_r,velx_r,vely_r,velz_r CCTK_REAL press_r,w_lorentz_r CCTK_REAL dens_l,tau_l,sx_l,sy_l,sz_l,rho_l,eps_l,velx_l,vely_l,velz_l CCTK_REAL press_l,w_lorentz_l CCTK_REAL dens_a,tau_a,sx_a,sy_a,sz_a,rho_a,eps_a,velx_a,vely_a,velz_a CCTK_REAL press_a,w_lorentz_a CCTK_REAL uxx,uxy,uxz,uyy,uyz,uzz,det,two CCTK_REAL gxxh,gxyh,gxzh,gyyh,gyzh,gzzh CCTK_REAL dens_vis,sx_vis,sy_vis,sz_vis,tau_vis CCTK_REAL dtx,dty,dtz,alpm,alpp,avg_alp,dxinv,dyinv,dzinv CCTK_REAL deps_by_drho,dpress_by_drho CCTK_REAL dx,dy,dz logical mahc_eos_ideal,mahc_eos_polytrope,mahc_eos_tables logical some_hydro_dust,some_hydro_mahc #ifdef CACTUSEOS_EOS_BASE #include "EOS_Base.inc" #endif some_hydro_dust = CCTK_EQUALS(some_hydro,'MahcHydroDust') some_hydro_mahc = CCTK_EQUALS(some_hydro,'MahcHydro') mahc_eos_ideal = CCTK_EQUALS(mahc_eos,'ideal_fluid') mahc_eos_polytrope = CCTK_EQUALS(mahc_eos,'polytrope') mahc_eos_tables = CCTK_EQUALS(mahc_eos,'tables') nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) c should we do TVD reconstuction? tvd = tvd_reconstruct c should we use eno? c if(CONTAINS('hydromethod','eno')) then c eno = .TRUE. c else c eno = .FALSE. c endif eno = .FALSE. dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) dtx = evolve_dt/dx dty = evolve_dt/dy dtz = evolve_dt/dz dxinv = 1.0d0/dx dyinv = 1.0d0/dy dzinv = 1.0d0/dz two = 2.0d0 c must get the slopes for the eno method c if (eno) then c if(order .eq. 2) then c call eno_slope_x(CCTK_FARGUMENTS) c else c eno_rho_slope = 0.0d0 c eno_eps_slope = 0.0d0 c eno_velx_slope = 0.0d0 c eno_vely_slope = 0.0d0 c eno_velz_slope = 0.0d0 c endif c endif c calculate numerical fluxes do k=1+mahc_stencil_size,nz-mahc_stencil_size do j=1+mahc_stencil_size,ny-mahc_stencil_size do i=mahc_stencil_size,nx-mahc_stencil_size c use averaged metric gxxh = 0.5d0 * (gxx(i,j,k)+gxx(i+1,j,k)) gxyh = 0.5d0 * (gxy(i,j,k)+gxy(i+1,j,k)) gxzh = 0.5d0 * (gxz(i,j,k)+gxz(i+1,j,k)) gyyh = 0.5d0 * (gyy(i,j,k)+gyy(i+1,j,k)) gyzh = 0.5d0 * (gyz(i,j,k)+gyz(i+1,j,k)) gzzh = 0.5d0 * (gzz(i,j,k)+gzz(i+1,j,k)) det= -(gxzh**2*gyyh) + two*gxyh*gxzh*gyzh - gxxh*gyzh**2 & - gxyh**2*gzzh + gxxh*gyyh*gzzh uxx=(-gyzh**2 + gyyh*gzzh)/det uxy=(gxzh*gyzh - gxyh*gzzh)/det uyy=(-gxzh**2 + gxxh*gzzh)/det uxz=(-gxzh*gyyh + gxyh*gyzh)/det uyz=(gxyh*gxzh - gxxh*gyzh)/det uzz=(-gxyh**2 + gxxh*gyyh)/det c calculate averaged shift, if we have it if(shift_state .eq. SHIFT_INACTIVE) then betaxh = 0.0d0 betayh = 0.0d0 betazh = 0.0d0 else betaxh = 0.5d0 * (betax(i,j,k)+betax(i+1,j,k)) betayh = 0.5d0 * (betay(i,j,k)+betay(i+1,j,k)) betazh = 0.5d0 * (betaz(i,j,k)+betaz(i+1,j,k)) endif #ifdef RECON_EVOL if(stencil_size .eq. 1) then write(*,*) 'flux_x_update: evolved reconstruction does not' write(*,*) 'work with stencil_size = 1!' STOP endif c get slopes and left and right states call slopelim(slopel,sloper,order,dens(i-1,j,k),dens(i,j,k), . dens(i+1,j,k),dens(i+2,j,k),tvd) dens_l = dens(i,j,k) + 0.5d0 * slopel*dxinv dens_r = dens(i+1,j,k) - 0.5d0 * sloper*dxinv call slopelim(slopel,sloper,order,tau(i-1,j,k),tau(i,j,k), . tau(i+1,j,k),tau(i+2,j,k),tvd) tau_l = tau(i,j,k) + 0.5d0 * slopel*dxinv tau_r = tau(i+1,j,k) - 0.5d0 * sloper*dxinv call slopelim(slopel,sloper,order,sx(i-1,j,k),sx(i,j,k), . sx(i+1,j,k),sx(i+2,j,k),tvd) sx_l = sx(i,j,k) + 0.5d0 * slopel*dxinv sx_r = sx(i+1,j,k) - 0.5d0 * sloper*dxinv call slopelim(slopel,sloper,order,sy(i-1,j,k),sy(i,j,k), . sy(i+1,j,k),sy(i+2,j,k),tvd) sy_l = sy(i,j,k) + 0.5d0 * slopel*dxinv sy_r = sy(i+1,j,k) - 0.5d0 * sloper*dxinv call slopelim(slopel,sloper,order,sz(i-1,j,k),sz(i,j,k), . sz(i+1,j,k),sz(i+2,j,k),tvd) sz_l = sz(i,j,k) + 0.5d0 * slopel*dxinv sz_r = sz(i+1,j,k) - 0.5d0 * sloper*dxinv c average evolved variables dens_a = 0.5d0 * (dens_l+dens_r) tau_a = 0.5d0 * (tau_l+tau_r) sx_a = 0.5d0 * (sx_l+sx_r) sy_a = 0.5d0 * (sy_l+sy_r) sz_a = 0.5d0 * (sz_l+sz_r) c get primative variables call hyp_get_prim_var_pt(dens_l,sx_l,sy_l,sz_l,tau_l,rho_l, . velx_l,vely_l,velz_l,eps_l,press_l,w_lorentz_l, . mahc_eos_gamma,uxx,uxy,uxz,uyy,uyz,uzz,det) call hyp_get_prim_var_pt(dens_r,sx_r,sy_r,sz_r,tau_r,rho_r, . velx_r,vely_r,velz_r,eps_r,press_r,w_lorentz_r, . mahc_eos_gamma,uxx,uxy,uxz,uyy,uyz,uzz,det) call hyp_get_prim_var_pt(dens_a,sx_a,sy_a,sz_a,tau_a,rho_a, . velx_a,vely_a,velz_a,eps_a,press_a,w_lorentz_a, . mahc_eos_gamma,uxx,uxy,uxz,uyy,uyz,uzz,det) #endif #ifdef RECON_PRIM c get slopes and left and right states if(eno) then c rho_l = rho(i ,j,k) + 0.5d0 * eno_rho_slope(i ,j,k) c rho_r = rho(i+1,j,k) - 0.5d0 * eno_rho_slope(i+1,j,k) c velx_l = velx(i ,j,k) + 0.5d0 * eno_velx_slope(i ,j,k) c velx_r = velx(i+1,j,k) - 0.5d0 * eno_velx_slope(i+1,j,k) c vely_l = vely(i ,j,k) + 0.5d0 * eno_vely_slope(i ,j,k) c vely_r = vely(i+1,j,k) - 0.5d0 * eno_vely_slope(i+1,j,k) c velz_l = velz(i ,j,k) + 0.5d0 * eno_velz_slope(i ,j,k) c velz_r = velz(i+1,j,k) - 0.5d0 * eno_velz_slope(i+1,j,k) else if(order .eq. 2) then call slopelim(slopel,sloper,order,rho(i-1,j,k),rho(i,j,k), . rho(i+1,j,k),rho(i+2,j,k),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif rho_l = rho(i,j,k) + 0.5d0 * slopel rho_r = rho(i+1,j,k) - 0.5d0 * sloper if(order .eq. 2) then call slopelim(slopel,sloper,order,velx(i-1,j,k),velx(i,j,k), . velx(i+1,j,k),velx(i+2,j,k),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif velx_l = velx(i,j,k) + 0.5d0 * slopel velx_r = velx(i+1,j,k) - 0.5d0 * sloper if(order .eq. 2) then call slopelim(slopel,sloper,order,vely(i-1,j,k),vely(i,j,k), . vely(i+1,j,k),vely(i+2,j,k),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif vely_l = vely(i,j,k) + 0.5d0 * slopel vely_r = vely(i+1,j,k) - 0.5d0 * sloper if(order .eq. 2) then call slopelim(slopel,sloper,order,velz(i-1,j,k),velz(i,j,k), . velz(i+1,j,k),velz(i+2,j,k),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif velz_l = velz(i,j,k) + 0.5d0 * slopel velz_r = velz(i+1,j,k) - 0.5d0 * sloper endif rho_a = 0.5d0 * (rho_l + rho_r) velx_a = 0.5d0 * (velx_l + velx_r) vely_a = 0.5d0 * (vely_l + vely_r) velz_a = 0.5d0 * (velz_l + velz_r) w_lorentz_l = 1.0d0 / sqrt(1.0d0 - ( . gxxh*velx_l**2 + gyyh*vely_l**2 + gzzh*velz_l**2 + . 2.0d0*gxyh*velx_l*vely_l + 2.0d0*gxzh*velx_l*velz_l + . 2.0d0*gyzh*vely_l*velz_l)) w_lorentz_r = 1.0d0 / sqrt(1.0d0 - ( . gxxh*velx_r**2 + gyyh*vely_r**2 + gzzh*velz_r**2 + . 2.0d0*gxyh*velx_r*vely_r + 2.0d0*gxzh*velx_r*velz_r + . 2.0d0*gyzh*vely_r*velz_r)) w_lorentz_a = 1.0d0 / sqrt(1.0d0 - ( . gxxh*velx_a**2 + gyyh*vely_a**2 + gzzh*velz_a**2 + . 2.0d0*gxyh*velx_a*vely_a + 2.0d0*gxzh*velx_a*velz_a + . 2.0d0*gyzh*vely_a*velz_a)) if(mahc_eos_ideal.or.mahc_eos_polytrope) then if(eno) then c eps_l = eps(i ,j,k) + 0.5d0 * eno_eps_slope(i ,j,k) c eps_r = eps(i+1,j,k) - 0.5d0 * eno_eps_slope(i+1,j,k) else if(order .eq. 2) then call slopelim(slopel,sloper,order,eps(i-1,j,k),eps(i,j,k), . eps(i+1,j,k),eps(i+2,j,k),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif eps_l = eps(i,j,k) + 0.5d0 * slopel eps_r = eps(i+1,j,k) - 0.5d0 * sloper endif press_l = (mahc_eos_gamma - 1.0d0) * rho_l * eps_l press_r = (mahc_eos_gamma - 1.0d0) * rho_r * eps_r if(mahc_eos_ideal) then eps_a = 0.5d0 * (eps_l + eps_r) press_a = (mahc_eos_gamma - 1.0d0) * rho_a * eps_a else eps_a = (mahc_eos_k * rho_a**(mahc_eos_gamma-1.0d0))/ . (mahc_eos_gamma-1.0d0) press_a = mahc_eos_k * rho_a**(mahc_eos_gamma) endif elseif (mahc_eos_tables) then #ifdef CACTUSEOS_EOS_BASE if(eno) then c eps_l = eps(i ,j,k) + 0.5d0 * eno_eps_slope(i ,j,k) c eps_r = eps(i+1,j,k) - 0.5d0 * eno_eps_slope(i+1,j,k) else if(order .eq. 2) then call slopelim(slopel,sloper,order,eps(i-1,j,k),eps(i,j,k), . eps(i+1,j,k),eps(i+2,j,k),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif eps_l = eps(i,j,k) + 0.5d0 * slopel eps_r = eps(i+1,j,k) - 0.5d0 * sloper endif press_l = EOS_Pressure(mahc_eos_handle,rho_l,eps_l) press_r = EOS_Pressure(mahc_eos_handle,rho_r,eps_r) eps_a = 0.5d0 * (eps_l + eps_r) press_a = EOS_Pressure(mahc_eos_handle,rho_a,eps_a) #else write(*,*) 'flux_x_update: if you want to use tables,' write(*,*) ' you must compile with thorn EOS_Base.' STOP #endif c old 1-d eos interface #ifdef THORN_EOS_GOTTA_REALLY_FIX_THIS_SOON call eos_tc_from_rho(rho_l,eps_l,press_l) call eos_tc_from_rho(rho_r,eps_r,press_r) call prim_var_eos_tc_from_rho(rho_a,eps_a,press_a, . deps_by_drho,dpress_by_drho) #endif endif dens_l = sqrt(det) * rho_l * w_lorentz_l sx_l = sqrt(det) * (rho_l * (1.0d0 + eps_l) + press_l) * . w_lorentz_l**2 * (gxxh*velx_l + gxyh*vely_l + gxzh*velz_l) sy_l = sqrt(det) * (rho_l * (1.0d0 + eps_l) + press_l) * . w_lorentz_l**2 * (gxyh*velx_l + gyyh*vely_l + gyzh*velz_l) sz_l = sqrt(det) * (rho_l * (1.0d0 + eps_l) + press_l) * . w_lorentz_l**2 * (gxzh*velx_l + gyzh*vely_l + gzzh*velz_l) tau_l = sqrt(det) * ((rho_l * (1.0d0 + eps_l) + press_l)* . w_lorentz_l**2 - press_l) - dens_l dens_r = sqrt(det) * rho_r * w_lorentz_r sx_r = sqrt(det) * (rho_r * (1.0d0 + eps_r) + press_r) * . w_lorentz_r**2 * (gxxh*velx_r + gxyh*vely_r + gxzh*velz_r) sy_r = sqrt(det) * (rho_r * (1.0d0 + eps_r) + press_r) * . w_lorentz_r**2 * (gxyh*velx_r + gyyh*vely_r + gyzh*velz_r) sz_r = sqrt(det) * (rho_r * (1.0d0 + eps_r) + press_r) * . w_lorentz_r**2 * (gxzh*velx_r + gyzh*vely_r + gzzh*velz_r) tau_r = sqrt(det) * ((rho_r * (1.0d0 + eps_r) + press_r)* . w_lorentz_r**2 - press_r) - dens_r dens_a = sqrt(det) * rho_a * w_lorentz_a sx_a = sqrt(det) * (rho_a * (1.0d0 + eps_a) + press_a) * . w_lorentz_a**2 * (gxxh*velx_a + gxyh*vely_a + gxzh*velz_a) sy_a = sqrt(det) * (rho_a * (1.0d0 + eps_a) + press_a) * . w_lorentz_a**2 * (gxyh*velx_a + gyyh*vely_a + gyzh*velz_a) sz_a = sqrt(det) * (rho_a * (1.0d0 + eps_a) + press_a) * . w_lorentz_a**2 * (gxzh*velx_a + gyzh*vely_a + gzzh*velz_a) tau_a = sqrt(det) * ((rho_a * (1.0d0 + eps_a) + press_a)* . w_lorentz_a**2 - press_a) - dens_a #endif c now, calculate the numerical flux avg_alp = 0.5d0 * (alp(i+1,j,k) + alp(i,j,k)) num_dens_flux(i,j,k) = . 0.5d0 * (dens_l*(velx_l - betaxh/avg_alp) + . dens_r*(velx_r - betaxh/avg_alp)) num_tau_flux(i,j,k) = . 0.5d0 * ( tau_l*(velx_l - betaxh/avg_alp) + . sqrt(det)*press_l*velx_l + . tau_r*(velx_r - betaxh/avg_alp) + . sqrt(det)*press_r*velx_r) num_sx_flux(i,j,k) = . 0.5d0 * ( sx_l*(velx_l - betaxh/avg_alp)+sqrt(det)*press_l + . sx_r*(velx_r - betaxh/avg_alp)+sqrt(det)*press_r) num_sy_flux(i,j,k) = . 0.5d0 * ( sy_l*(velx_l - betaxh/avg_alp) + . sy_r*(velx_r - betaxh/avg_alp)) num_sz_flux(i,j,k) = . 0.5d0 * ( sz_l*(velx_l - betaxh/avg_alp) + . sz_r*(velx_r - betaxh/avg_alp)) c now, calculate the viscosity part of the numerical flux #ifndef CACTUSEOS_EOS_BASE if(mahc_eos_ideal) then dpdrho = (mahc_eos_gamma - 1.0d0) * eps_a dpdeps = (mahc_eos_gamma - 1.0d0) * rho_a fmi = 0 else dpdrho = mahc_eos_k * mahc_eos_gamma * rho_a**(mahc_eos_gamma-1.0d0) dpdeps = 0.0d0 fmi = 1 endif #else if(mahc_eos_ideal) then dpdrho = (mahc_eos_gamma - 1.0d0) * eps_a dpdeps = (mahc_eos_gamma - 1.0d0) * rho_a fmi = 0 elseif(mahc_eos_polytrope) then dpdrho = mahc_eos_k * mahc_eos_gamma * rho_a**(mahc_eos_gamma-1.0d0) dpdeps = 0.0d0 fmi = 1 elseif(mahc_eos_tables) then dpdrho = EOS_DPressByDRho(mahc_eos_handle,rho_a,eps_a) dpdeps = EOS_DPressByDEps(mahc_eos_handle,rho_a,eps_a) fmi = 0 endif c elseif(mahc_eos_tables) then c dpdrho = dpress_by_drho c dpdeps = 0.0d0 c fmi = 1 c endif #endif call vis_contrib(dens_vis,sx_vis,sy_vis,sz_vis,tau_vis, . dens_a,sx_a,sy_a,sz_a,tau_a,rho_a,velx_a,vely_a,velz_a, . eps_a,press_a,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, . uxx,uxy,uxz,uyy,uyz,uzz,det,dpdrho,dpdeps,avg_alp, . dens_r,sx_r,sy_r,sz_r,tau_r,dens_l,sx_l,sy_l,sz_l,tau_l, . betaxh,betayh,betazh,fmi,some_hydro_dust,some_hydro_mahc) c add in viscosity part num_dens_flux(i,j,k) = num_dens_flux(i,j,k) + dens_vis num_sx_flux(i,j,k) = num_sx_flux(i,j,k) + sx_vis num_sy_flux(i,j,k) = num_sy_flux(i,j,k) + sy_vis num_sz_flux(i,j,k) = num_sz_flux(i,j,k) + sz_vis num_tau_flux(i,j,k) = num_tau_flux(i,j,k) + tau_vis enddo enddo enddo c calculate the updates! do k=1+mahc_stencil_size,nz-mahc_stencil_size do j=1+mahc_stencil_size,ny-mahc_stencil_size do i=1+mahc_stencil_size,nx-mahc_stencil_size alpm = 0.5d0 * (alp(i-1,j,k)+alp(i,j,k)) alpp = 0.5d0 * (alp(i+1,j,k)+alp(i,j,k)) del_dens(i,j,k) = del_dens(i,j,k) -dtx * . (alpp*num_dens_flux(i,j,k) - alpm*num_dens_flux(i-1,j,k)) del_sx(i,j,k) = del_sx(i,j,k) - dtx * . (alpp*num_sx_flux(i,j,k) - alpm*num_sx_flux(i-1,j,k)) del_sy(i,j,k) = del_sy(i,j,k) - dtx * . (alpp*num_sy_flux(i,j,k) - alpm*num_sy_flux(i-1,j,k)) del_sz(i,j,k) = del_sz(i,j,k) - dtx * . (alpp*num_sz_flux(i,j,k) - alpm*num_sz_flux(i-1,j,k)) del_tau(i,j,k) = del_tau(i,j,k) -dtx * . (alpp*num_tau_flux(i,j,k) - alpm*num_tau_flux(i-1,j,k)) enddo enddo enddo return end