c Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations c Copyright (C) 2001 Mark Miller /*@@ @file roe_flux_z_update.F @date December 1999 @author Mark Miller @desc Performs the Z 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_z_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_z_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 dtz,alpm,alpp,avg_alp,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) dtz = evolve_dt/dz 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_z(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=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 c use averaged metric gxxh = 0.5d0 * (gxx(i,j,k)+gxx(i,j,k+1)) gxyh = 0.5d0 * (gxy(i,j,k)+gxy(i,j,k+1)) gxzh = 0.5d0 * (gxz(i,j,k)+gxz(i,j,k+1)) gyyh = 0.5d0 * (gyy(i,j,k)+gyy(i,j,k+1)) gyzh = 0.5d0 * (gyz(i,j,k)+gyz(i,j,k+1)) gzzh = 0.5d0 * (gzz(i,j,k)+gzz(i,j,k+1)) 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,j,k+1)) betayh = 0.5d0 * (betay(i,j,k)+betay(i,j,k+1)) betazh = 0.5d0 * (betaz(i,j,k)+betaz(i,j,k+1)) 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,j,k-1),dens(i,j,k), . dens(i,j,k+1),dens(i,j,k+2),tvd) dens_l = dens(i,j,k) + 0.5d0 * slopel*dzinv dens_r = dens(i,j,k+1) - 0.5d0 * sloper*dzinv call slopelim(slopel,sloper,order,tau(i,j,k-1),tau(i,j,k), . tau(i,j,k+1),tau(i,j,k+2),tvd) tau_l = tau(i,j,k) + 0.5d0 * slopel*dzinv tau_r = tau(i,j,k+1) - 0.5d0 * sloper*dzinv call slopelim(slopel,sloper,order,sx(i,j,k-1),sx(i,j,k), . sx(i,j,k+1),sx(i,j,k+2),tvd) sx_l = sx(i,j,k) + 0.5d0 * slopel*dzinv sx_r = sx(i,j,k+1) - 0.5d0 * sloper*dzinv call slopelim(slopel,sloper,order,sy(i,j,k-1),sy(i,j,k), . sy(i,j,k+1),sy(i,j,k+2),tvd) sy_l = sy(i,j,k) + 0.5d0 * slopel*dzinv sy_r = sy(i,j,k+1) - 0.5d0 * sloper*dzinv call slopelim(slopel,sloper,order,sz(i,j,k-1),sz(i,j,k), . sz(i,j,k+1),sz(i,j,k+2),tvd) sz_l = sz(i,j,k) + 0.5d0 * slopel*dzinv sz_r = sz(i,j,k+1) - 0.5d0 * sloper*dzinv 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,j,k+1) - 0.5d0 * eno_rho_slope(i,j,k+1) c velx_l = velx(i,j,k ) + 0.5d0 * eno_velx_slope(i,j,k ) c velx_r = velx(i,j,k+1) - 0.5d0 * eno_velx_slope(i,j,k+1) c vely_l = vely(i,j,k ) + 0.5d0 * eno_vely_slope(i,j,k ) c vely_r = vely(i,j,k+1) - 0.5d0 * eno_vely_slope(i,j,k+1) c velz_l = velz(i,j,k ) + 0.5d0 * eno_velz_slope(i,j,k ) c velz_r = velz(i,j,k+1) - 0.5d0 * eno_velz_slope(i,j,k+1) else if(order .eq. 2) then call slopelim(slopel,sloper,order,rho(i,j,k-1),rho(i,j,k), . rho(i,j,k+1),rho(i,j,k+2),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif rho_l = rho(i,j,k) + 0.5d0 * slopel rho_r = rho(i,j,k+1) - 0.5d0 * sloper if(order .eq. 2) then call slopelim(slopel,sloper,order,velx(i,j,k-1),velx(i,j,k), . velx(i,j,k+1),velx(i,j,k+2),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif velx_l = velx(i,j,k) + 0.5d0 * slopel velx_r = velx(i,j,k+1) - 0.5d0 * sloper if(order .eq. 2) then call slopelim(slopel,sloper,order,vely(i,j,k-1),vely(i,j,k), . vely(i,j,k+1),vely(i,j,k+2),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif vely_l = vely(i,j,k) + 0.5d0 * slopel vely_r = vely(i,j,k+1) - 0.5d0 * sloper if(order .eq. 2) then call slopelim(slopel,sloper,order,velz(i,j,k-1),velz(i,j,k), . velz(i,j,k+1),velz(i,j,k+2),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif velz_l = velz(i,j,k) + 0.5d0 * slopel velz_r = velz(i,j,k+1) - 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,j,k+1) - 0.5d0 * eno_eps_slope(i,j,k+1) else if(order .eq. 2) then call slopelim(slopel,sloper,order,eps(i,j,k-1),eps(i,j,k), . eps(i,j,k+1),eps(i,j,k+2),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif eps_l = eps(i,j,k) + 0.5d0 * slopel eps_r = eps(i,j,k+1) - 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,j,k+1) - 0.5d0 * eno_eps_slope(i,j,k+1) else if(order .eq. 2) then call slopelim(slopel,sloper,order,eps(i,j,k-1),eps(i,j,k), . eps(i,j,k+1),eps(i,j,k+2),tvd) else slopel = 0.0d0 sloper = 0.0d0 endif eps_l = eps(i,j,k) + 0.5d0 * slopel eps_r = eps(i,j,k+1) - 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_z_update: if you want to use tables,' write(*,*) ' you must compile with thorn EOS_Base.' STOP #endif c old 1-d eos interface c#ifdef THORN_EOS_GOTTA_REALLY_FIX_THIS_SOON c call eos_tc_from_rho(rho_l,eps_l,press_l) c call eos_tc_from_rho(rho_r,eps_r,press_r) c call prim_var_eos_tc_from_rho(rho_a,eps_a,press_a, c . deps_by_drho,dpress_by_drho) 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,j,k+1) + alp(i,j,k)) num_dens_flux(i,j,k) = . 0.5d0 * (dens_l*(velz_l - betazh/avg_alp) + . dens_r*(velz_r - betazh/avg_alp)) num_tau_flux(i,j,k) = . 0.5d0 * ( tau_l*(velz_l - betazh/avg_alp) + . sqrt(det)*press_l*velz_l + . tau_r*(velz_r - betazh/avg_alp) + . sqrt(det)*press_r*velz_r) num_sx_flux(i,j,k) = . 0.5d0 * ( sx_l*(velz_l - betazh/avg_alp) + . sx_r*(velz_r - betazh/avg_alp)) num_sy_flux(i,j,k) = . 0.5d0 * ( sy_l*(velz_l - betazh/avg_alp) + . sy_r*(velz_r - betazh/avg_alp)) num_sz_flux(i,j,k) = . 0.5d0 * ( sz_l*(velz_l - betazh/avg_alp)+sqrt(det)*press_l + . sz_r*(velz_r - betazh/avg_alp)+sqrt(det)*press_r) 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 c elseif(mahc_eos_tables) then c dpdrho = dpress_by_drho c dpdeps = 0.0d0 c fmi = 1 endif #endif call vis_contrib(dens_vis,sz_vis,sx_vis,sy_vis,tau_vis, . dens_a,sz_a,sx_a,sy_a,tau_a,rho_a,velz_a,velx_a,vely_a, . eps_a,press_a,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh, . uzz,uxz,uyz,uxx,uxy,uyy,det,dpdrho,dpdeps,avg_alp, . dens_r,sz_r,sx_r,sy_r,tau_r,dens_l,sz_l,sx_l,sy_l,tau_l, . betazh,betaxh,betayh,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,j,k-1)+alp(i,j,k)) alpp = 0.5d0 * (alp(i,j,k+1)+alp(i,j,k)) del_dens(i,j,k) = del_dens(i,j,k) -dtz * . (alpp*num_dens_flux(i,j,k) - alpm*num_dens_flux(i,j,k-1)) del_sx(i,j,k) = del_sx(i,j,k) - dtz * . (alpp*num_sx_flux(i,j,k) - alpm*num_sx_flux(i,j,k-1)) del_sy(i,j,k) = del_sy(i,j,k) - dtz * . (alpp*num_sy_flux(i,j,k) - alpm*num_sy_flux(i,j,k-1)) del_sz(i,j,k) = del_sz(i,j,k) - dtz * . (alpp*num_sz_flux(i,j,k) - alpm*num_sz_flux(i,j,k-1)) del_tau(i,j,k) = del_tau(i,j,k) -dtz * . (alpp*num_tau_flux(i,j,k) - alpm*num_tau_flux(i,j,k-1)) enddo enddo enddo return end