diff options
| author | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-01-14 14:23:37 +0000 |
|---|---|---|
| committer | rhaas <rhaas@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2013-01-14 14:23:37 +0000 |
| commit | 5199c61af4721444a26c4e7c1c30d3b8fdedf54c (patch) | |
| tree | 9428f573d53e23f6b2970cbe664ba96cfb2eb2ec /src | |
| parent | 48b053d41b2aa9ce720d443d288eef535fb9651d (diff) | |
GRHydro: Add basic vector potential support
Basic cell-centered, algebraic gauge vector potential
method with place-holders for lorenz gauge.
Initial Avec constrained to poloidal at the moment.
From: Tanja Bode <tanja.bode@physics.gatech.edu>
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@456 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
Diffstat (limited to 'src')
| -rw-r--r-- | src/GRHydro_BvecfromAvec.F90 | 147 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimAM.F90 | 1365 | ||||
| -rw-r--r-- | src/GRHydro_Con2PrimM.F90 | 2 | ||||
| -rw-r--r-- | src/GRHydro_HLLE_AM.F90 | 943 | ||||
| -rw-r--r-- | src/GRHydro_InterfacesAM.h | 95 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheck.F90 | 10 | ||||
| -rw-r--r-- | src/GRHydro_Prim2ConAM.F90 | 762 | ||||
| -rw-r--r-- | src/GRHydro_RegisterVars.cc | 6 | ||||
| -rw-r--r-- | src/GRHydro_RiemannSolveAM.F90 | 144 | ||||
| -rw-r--r-- | src/GRHydro_SourceAM.F90 | 679 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMaskM.F90 | 475 | ||||
| -rw-r--r-- | src/make.code.defn | 6 |
12 files changed, 4632 insertions, 2 deletions
diff --git a/src/GRHydro_BvecfromAvec.F90 b/src/GRHydro_BvecfromAvec.F90 new file mode 100644 index 0000000..ea0a58d --- /dev/null +++ b/src/GRHydro_BvecfromAvec.F90 @@ -0,0 +1,147 @@ + /*@@ + @file GRHydro_BvecfromAvec + @date Aug 31, 2010 + @author Tanja Bode + @desc + Calculate B^i (at cell center) from Avec (at edges) + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Avecx(i,j,k) Avec(i,j,k,1) +#define Avecy(i,j,k) Avec(i,j,k,2) +#define Avecz(i,j,k) Avec(i,j,k,3) + +! Second order f.d. for volume-centered quantity from volume-centered quantity +#define DIFF_X_2(q) (0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx) +#define DIFF_Y_2(q) (0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy) +#define DIFF_Z_2(q) (0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz) + +! Fourth order f.d. for volume-centered quantity from volume-centered quantity +#define DIFF_X_4(q) ((-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ + q(i-2,j,k)) / 12.d0 * idx) +#define DIFF_Y_4(q) ((-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ + q(i,j-2,k)) / 12.d0 * idy) +#define DIFF_Z_4(q) ((-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ + q(i,j,k-2)) / 12.d0 * idz) + + /*@@ + @routine GRHydro_BvecfromAvec + @date Oct 24 + @author Tanja Bode + @desc + Calculate B^i (at cell center) from Avec (at edges) + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_BvecfromAvec(CCTK_ARGUMENTS) + + implicit none + + CCTK_REAL :: Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y + CCTK_REAL :: det, isdet + CCTK_REAL :: dx, dy, dz, idx, idy, idz + CCTK_INT :: i, j, k, nx, ny, nz, GRHydro_UseGeneralCoordinates + + logical, allocatable, dimension (:,:,:) :: force_spatial_second_order + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + !!DECLARE_CCTK_FUNCTIONS + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + call CCTK_WARN(0,"Bvec from Avec only in Cartesian at the moment."); + end if + + 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) + idx = 1.d0/dx + idy = 1.d0/dy + idz = 1.d0/dz + + allocate (force_spatial_second_order(nx,ny,nz)) + force_spatial_second_order = .FALSE. + + if (spatial_order > 2) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = 1 + GRHydro_stencil, nz - GRHydro_stencil + do j = 1 + GRHydro_stencil, ny - GRHydro_stencil + do i = 1 + GRHydro_stencil, nx - GRHydro_stencil + if ((i < 3).or.(i > cctk_lsh(1) - 2).or. & + (j < 3).or.(j > cctk_lsh(2) - 2).or. & + (k < 3).or.(k > cctk_lsh(3) - 2) ) then + force_spatial_second_order(i,j,k) = .TRUE. + else if ( use_mask > 0 ) then + if (minval(emask(i-2:i+2,j-2:j+2,k-2:k+2)) < 0.75d0) then + force_spatial_second_order(i,j,k) = .TRUE. + end if + end if + end do + end do + end do + !$OMP END PARALLEL DO + end if + + + call CCTK_WARN(1,"Bvec from Avec start."); + !$OMP PARALLEL DO PRIVATE( det, isdet, local_spatial_order, & + !$OMP Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y) + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + + local_spatial_order = spatial_order + if (force_spatial_second_order(i,j,k)) then + local_spatial_order = 2 + end if + + if (local_spatial_order .eq. 2) then + Az_y = DIFF_Y_2(Avecz) + Ay_z = DIFF_Z_2(Avecy) + Ax_z = DIFF_Z_2(Avecx) + Az_x = DIFF_X_2(Avecz) + Ay_x = DIFF_X_2(Avecy) + Ax_y = DIFF_Y_2(Avecx) + else + Az_y = DIFF_Y_4(Avecz) + Ay_z = DIFF_Z_4(Avecy) + Ax_z = DIFF_Z_4(Avecx) + Az_x = DIFF_X_4(Avecz) + Ay_x = DIFF_X_4(Avecy) + Ax_y = DIFF_Y_4(Avecx) + end if + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + isdet = 1.d0/sqrt(det) + + Bvecx(i,j,k) = isdet*( Az_y - Ay_z ) + Bvecy(i,j,k) = isdet*( Ax_z - Az_x ) + Bvecz(i,j,k) = isdet*( Ay_x - Ax_y ) + + end do + end do + end do + !$OMP END PARALLEL DO + call CCTK_WARN(1,"Bvec from Avec end."); + + +end subroutine GRHydro_BvecfromAvec + diff --git a/src/GRHydro_Con2PrimAM.F90 b/src/GRHydro_Con2PrimAM.F90 new file mode 100644 index 0000000..50d701c --- /dev/null +++ b/src/GRHydro_Con2PrimAM.F90 @@ -0,0 +1,1365 @@ +/*@@ + @file GRHydro_Con2PrimAM.F90 + @date Sep 3, 2010 + @author Scott Noble, Joshua Faber, Bruno Mundim, Tanja Bode + @desc + The routines for converting conservative to primitive variables. + Vector-potential MHD. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "SpaceMask.h" +#include "GRHydro_InterfacesAM.h" +#include "GRHydro_Macros.h" + +#define ITER_TOL (1.0e-8) +#define MAXITER (50) + + /*@@ + @routine Conservative2PrimitiveAM + @date Sep 3, 2010 + @author Scott Noble, Joshua Faber, Bruno Mundim, Ian Hawke + @desc + Wrapper routine that converts from conserved to primitive variables + at every grid cell centre. + @enddesc + @calls + @calledby + @history + Trimmed and altered from the GR3D routines, original author Mark Miller. + 2007?: Bruno excluded the points in the atmosphere and excision region from the computation. + Aug. 2008: Luca added a check on whether a failure at a given point may be disregarded, + because that point will then be restricted from a finer level. This should be completely + safe only if *regridding happens at times when all levels are evolved.* + Feb. 2009: The above procedure proved to be wrong, so Luca implemented another one. + When a failure occurs, it is temporarily ignored, except for storing the location of where + it occured in a mask. At the end, after a Carpet restriction, the mask is checked and if + it still contains failures, the run is aborted with an error message. Only used routines + have been updated to use this procedure. + @endhistory + +@@*/ + +subroutine Conservative2PrimitiveAM(CCTK_ARGUMENTS) + + use Con2PrimM_fortran_interfaces + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: i, j, k, itracer, nx, ny, nz + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, sdet, pmin(1), epsmin(1) + CCTK_REAL :: oob, b2, d2, s2, bscon, bxhat, byhat, bzhat, bhatscon + CCTK_REAL :: Wm, Wm0, Wm_resid, Wmold + CCTK_REAL :: s2m, s2m0, s2m_resid, s2mold, s2max, taum + CCTK_INT :: niter + CCTK_INT :: epsnegative + character(len=100) warnline + + CCTK_REAL :: local_min_tracer, local_gam(1), local_pgam,local_K, sc + CCTK_REAL :: local_perc_ptol + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1),one(1)=1.0d0 +! end EOS Omni vars + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + + logical :: posdef + + CCTK_REAL :: g11c, g12c, g13c, g22c, g23c, g33c + CCTK_REAL :: tmp1 + + ! Save the primitive variables to temporary functions before calling the + ! con2prim pointwise routines: + CCTK_REAL :: rho_tmp, press_tmp, eps_tmp + CCTK_REAL :: velx_tmp, vely_tmp, velz_tmp, w_lorentz_tmp + CCTK_REAL :: Bvecx_tmp, Bvecy_tmp, Bvecz_tmp + CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp + + CCTK_REAL :: bdotv, magpress + + ! Assume 3-metric is positive definite. Check deep inside the horizon + ! if this is actually satisfied and if it is not then cast the metric + !as conformally flat only for con2prim inversion purposes. + posdef = .true. + + !! Get Bvec from Avec before rest, so we can reuse pointwise MHD recovery + call GRHydro_BvecfromAvec(CCTK_PASS_FTOF) + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define vel faulty_vel +#define Bvec faulty_Bvec + +! begin EOS Omni vars + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0;xeps(1)=0.0d0 +! end EOS Omni vars + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + if (use_min_tracer .ne. 0) then + local_min_tracer = min_tracer + else + local_min_tracer = 0d0 + end if + +! call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& +! GRHydro_rho_min,xeps,xtemp,xye,pmin,keyerr,anyerr) +! call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& +! GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) + ! this is a poly call + xrho(1)=GRHydro_rho_min + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,pmin,keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + one,one,xtemp,xye,local_gam,keyerr,anyerr) + local_gam = local_gam+1.d0 + + !call CCTK_WARN(1,"In Con2PrimAM") + + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & + !$OMP b2,xrho,xeps,xpress,xtemp,local_K,local_pgam,sc,keyerr,anyerr,keytemp, & + !$OMP local_perc_ptol,posdef,g11c,g12c,g13c,g22c,g23c,g33c,tmp1, & + !$OMP sdet,d2,s2,oob,bscon,bxhat,byhat,bzhat, & + !$OMP bhatscon,Wm,Wm0,Wm_resid,Wmold,s2m,s2m0,s2m_resid,s2mold,s2max, & + !$OMP taum,niter,rho_tmp,press_tmp,eps_tmp,velx_tmp,vely_tmp,velz_tmp, & + !$OMP Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, & + !$OMP w_lorentz_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,bdotv,magpress) + do k = 1, nz + do j = 1, ny + do i = 1, nx + + !do not compute if in atmosphere or in excised region + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & + (hydro_excision_mask(i,j,k) .ne. 0)) cycle + + epsnegative = 0 + + det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) + sdet = sqrt(det) + + + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& + g23(i,j,k),g33(i,j,k)) + +!!$ Tracers don't need an MHD treatment! + if (evolve_tracer .ne. 0) then + do itracer=1,number_of_tracers + call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), tracer(i,j,k,itracer), & + dens(i,j,k)) + + if (use_min_tracer .ne. 0) then + if (tracer(i,j,k,itracer) .le. local_min_tracer) then + tracer(i,j,k,itracer) = local_min_tracer + end if + end if + + enddo + + endif + + if(evolve_Y_e.ne.0) then + Y_e(i,j,k) = max(min(Y_e_con(i,j,k) / dens(i,j,k),GRHydro_Y_e_max),& + GRHydro_Y_e_min) + endif + + + b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ & + 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ & + g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3)) + + + if ( dens(i,j,k) .le. sdet*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + + !call CCTK_WARN(1,"Con2Prim: Resetting to atmosphere") + !write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + !call CCTK_WARN(1,warnline) + !write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k),& + ! temperature(i,j,k),y_e(i,j,k) + !call CCTK_WARN(1,warnline) + + + dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + rho(i,j,k) = GRHydro_rho_min + scon(i,j,k,:) = 0.d0 + vup(i,j,k,:) = 0.d0 + w_lorentz(i,j,k) = 1.d0 + + if(evolve_temper.ne.0) then + ! set hot atmosphere values + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k) + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr) + else + keytemp = 0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + endif + + !call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + ! rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + ! + !call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + ! rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + + ! w_lorentz=1, so the expression for tau reduces to: + + !!$ tau does need to take into account the existing B-field + !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] + tau(i,j,k) = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + + if(tau(i,j,k).le.sdet*b2*0.5d0)then + tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0 + endif + + cycle + + end if + + + if(evolve_temper.eq.0) then + + + if(sqrtdet_thr.gt.0d0 .and. sdet.ge.sqrtdet_thr) then + d2 = dens(i,j,k)**2 + s2 = uxx*scon(i,j,k,1)**2 + uyy*scon(i,j,k,2)**2 & + + uzz*scon(i,j,k,3)**2 & + + 2.0d0*uxy*scon(i,j,k,1)*scon(i,j,k,2) & + + 2.0d0*uxz*scon(i,j,k,1)*scon(i,j,k,3) & + + 2.0d0*uyz*scon(i,j,k,2)*scon(i,j,k,3) + oob = 1.0d0/sqrt(b2) + bxhat = oob*Bprim(i,j,k,1) + byhat = oob*Bprim(i,j,k,2) + bzhat = oob*Bprim(i,j,k,3) + bhatscon = bxhat*scon(i,j,k,1)+byhat*scon(i,j,k,2) & + +bzhat*scon(i,j,k,3) + bscon = Bprim(i,j,k,1)*scon(i,j,k,1) & + + Bprim(i,j,k,2)*scon(i,j,k,2) & + + Bprim(i,j,k,3)*scon(i,j,k,3) + ! Initial guesses for iterative procedure to find Wm: + Wm0 = sdet*sqrt(bhatscon**2+d2) + s2m0 = (Wm0**2*s2+bhatscon**2*(b2+2.0d0*Wm0)) & + / (Wm0+b2)**2 + Wm = sdet*sqrt(s2m0+d2) + s2m = (Wm**2*s2+bscon**2*(b2+2.0d0*Wm)) & + / (Wm+b2)**2 + s2m_resid = 1.0d60 + Wm_resid = 1.0d60 + niter = 0 + do while((s2m_resid.ge.ITER_TOL.and.Wm_resid.ge.ITER_TOL).and.& + niter.le.MAXITER) + Wmold = Wm + s2mold = s2m + Wm = sdet*sqrt(s2m+d2) + s2m = (Wm**2*s2+bscon**2*(b2+2.0d0*Wm)) & + / (Wm+b2)**2 + Wm_resid = abs(Wmold-Wm) + s2m_resid = abs(s2mold-s2m) + niter = niter + 1 + end do + !TODO: abort execution if niter .eq. MAXITER and warn user + taum = tau(i,j,k) - 0.5d0*sdet*b2 -0.5d0*(b2*s2-bscon**2)/ & + (sdet*(Wm+b2)**2) + s2max = taum*(taum+2.0d0*dens(i,j,k)) + if(taum.lt.GRHydro_tau_min)then + tau(i,j,k) = GRHydro_tau_min + 0.5d0*sdet*b2 + 0.5d0* & + (b2*s2-bscon**2)/(sdet*(Wm+b2)**2) + end if + if(s2.gt.s2max) then + scon(i,j,k,1) = scon(i,j,k,1)*sqrt(s2max/s2) + scon(i,j,k,2) = scon(i,j,k,2)*sqrt(s2max/s2) + scon(i,j,k,3) = scon(i,j,k,3)*sqrt(s2max/s2) + end if + endif + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + + !! We've already create B from A, construct Bcons + Bconsx_tmp = sdet*Bvecx_tmp + Bconsy_tmp = sdet*Bvecy_tmp + Bconsz_tmp = sdet*Bvecz_tmp + + keytemp = 0 + + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,xye(1), & + xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + if(evolve_entropy.ne.0) then + entropy(i,j,k) = entropycons(i,j,k)/dens(i,j,k)*rho(i,j,k) + endif + + if(GRHydro_C2P_failed(i,j,k).ne.0) then + xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress(1); + + xrho=10.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress(1)/local_K)/log(xrho(1)) + sc = local_K*dens(i,j,k) + + if(sdet.ge.sqrtdet_thr) then + GRHydro_C2P_failed(i,j,k) = 0 + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + Bconsx_tmp = sdet*Bvecx_tmp + Bconsy_tmp = sdet*Bvecy_tmp + Bconsz_tmp = sdet*Bvecz_tmp + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + if(GRHydro_C2P_failed(i,j,k).ne.0) then + GRHydro_C2P_failed(i,j,k) = 0 + call prim2conAM(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k), & + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),det, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), & + tau(i,j,k), & + rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3), & + eps(i,j,k),press(i,j,k), & + Bprim(i,j,k,1),Bprim(i,j,k,2),Bprim(i,j,k,3),w_lorentz(i,j,k)) + cycle + end if + end if + + bdotv=g11(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,1)+ & + g22(i,j,k)*Bprim(i,j,k,2)*vup(i,j,k,2)+ & + g33(i,j,k)*Bprim(i,j,k,3)*vup(i,j,k,3)+ & + 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,2)+ & + g13(i,j,k)*Bprim(i,j,k,1)*vup(i,j,k,3)+ & + g23(i,j,k)*Bprim(i,j,k,2)*vup(i,j,k,3)) + + magpress = 0.5d0*(b2/w_lorentz(i,j,k)**2+bdotv**2) + + if(rho(i,j,k)*eps(i,j,k)*max_magnetic_to_gas_pressure_ratio.le.magpress) then + GRHydro_C2P_failed(i,j,k) = 0 + + if(evolve_entropy.ne.0) then + local_K = entropycons(i,j,k)/dens(i,j,k) + xrho=10.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress(1)/local_K)/log(xrho(1)) + sc = local_K*dens(i,j,k) + end if + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + Bconsx_tmp = sdet*Bvecx_tmp + Bconsy_tmp = sdet*Bvecy_tmp + Bconsz_tmp = sdet*Bvecz_tmp + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + rho(i,j,k) = rho_tmp + press(i,j,k) = press_tmp + eps(i,j,k) = eps_tmp + vup(i,j,k,1) = velx_tmp + vup(i,j,k,2) = vely_tmp + vup(i,j,k,3) = velz_tmp + w_lorentz(i,j,k) = w_lorentz_tmp + Bprim(i,j,k,1) = Bvecx_tmp + Bprim(i,j,k,2) = Bvecy_tmp + Bprim(i,j,k,3) = Bvecz_tmp + cycle + end if + end if + + else ! if(evolve_temper.eq.0) then + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + Bconsx_tmp = sdet*Bvecx_tmp + Bconsy_tmp = sdet*Bvecy_tmp + Bconsz_tmp = sdet*Bvecz_tmp + + keytemp = 0 + + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & + GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, Y_e(i,j,k), & + temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + if(GRHydro_C2P_failed(i,j,k).ne.0) then + ! this means c2p did not converge. + ! In this case, we attempt to call c2p with a reduced + ! accuracy requirement; if it fails again, we abort + GRHydro_C2P_failed(i,j,k) = 0 + local_perc_ptol = GRHydro_eos_rf_prec*100.0d0 + ! Use the previous primitive values as initial guesses + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + Bconsx_tmp = sdet*Bvecx_tmp + Bconsy_tmp = sdet*Bvecy_tmp + Bconsz_tmp = sdet*Bvecz_tmp + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, & + local_perc_ptol, local_gam(1), dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, Y_e(i,j,k), & + temperature(i,j,k),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& + eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& + w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& + g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + if(GRHydro_C2P_failed(i,j,k).ne.0) then + !$OMP CRITICAL + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + call CCTK_WARN(1,"Convergence problem in c2p") + write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,"(3i5,1P10E15.6)") i,j,k,x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k),& + temperature(i,j,k),y_e(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif + !$OMP END CRITICAL + endif + endif + endif ! if(evolve_temper.eq.0) then + + + if (epsnegative .ne. 0) then + +#if 0 + ! cott 2010/03/30: + ! Set point to atmosphere, but continue evolution -- this is better than using + ! the poly EOS -- it will lead the code to crash if this happens inside a (neutron) star, + ! but will allow the job to continue if it happens in the atmosphere or in a + ! zone that contains garbage (i.e. boundary, buffer zones) + ! Ultimately, we want this fixed via a new carpet mask presently under development + ! GRHydro_C2P_failed(i,j,k) = 1 + + !$OMP CRITICAL + call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0! ') + write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel + call CCTK_WARN(GRHydro_NaN_verbose+2,warnline) + write(warnline,'(a20,3g16.7)') 'xyz location: ',& + x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(GRHydro_NaN_verbose+2,warnline) + write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) + call CCTK_WARN(GRHydro_NaN_verbose+2,warnline) + call CCTK_WARN(GRHydro_NaN_verbose+2,"Setting the point to atmosphere") + !$OMP END CRITICAL + + ! for safety, let's set the point to atmosphere + dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + rho(i,j,k) = GRHydro_rho_min + scon(i,j,k,:) = 0.d0 + vup(i,j,k,:) = 0.d0 + w_lorentz(i,j,k) = 1.d0 + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + + b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ & + 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ & + g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3)) + + + ! w_lorentz=1, so the expression for tau reduces to [see above]: + tau(i,j,k) = sqrt(det) * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) +#else + ! cott 2010/03/27: + ! Honestly, this should never happen. We need to flag the point where + ! this happened as having led to failing con2prim. + + !$OMP CRITICAL + call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype.') + !$OMP END CRITICAL + + xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress(1); + + xrho=10.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress(1)/local_K)/log(xrho(1)) + sc = local_K*dens(i,j,k) + + rho_tmp = rho(i,j,k) + press_tmp = press(i,j,k) + eps_tmp = eps(i,j,k) + velx_tmp = vup(i,j,k,1) + vely_tmp = vup(i,j,k,2) + velz_tmp = vup(i,j,k,3) + w_lorentz_tmp = w_lorentz(i,j,k) + Bvecx_tmp = Bprim(i,j,k,1) + Bvecy_tmp = Bprim(i,j,k,2) + Bvecz_tmp = Bprim(i,j,k,3) + Bconsx_tmp = sdet*Bvecx_tmp + Bconsy_tmp = sdet*Bvecy_tmp + Bconsz_tmp = sdet*Bvecz_tmp + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp,rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + +#endif + + end if ! if (epsnegative .ne. 0) then + + rho(i,j,k) = rho_tmp + press(i,j,k) = press_tmp + eps(i,j,k) = eps_tmp + vup(i,j,k,1) = velx_tmp + vup(i,j,k,2) = vely_tmp + vup(i,j,k,3) = velz_tmp + w_lorentz(i,j,k) = w_lorentz_tmp + Bprim(i,j,k,1) = Bvecx_tmp + Bprim(i,j,k,2) = Bvecy_tmp + Bprim(i,j,k,3) = Bvecz_tmp + + if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance)) then +! if ( rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) .or. GRHydro_C2P_failed(i,j,k) .ge. 1) then + + b2=g11(i,j,k)*Bprim(i,j,k,1)**2+g22(i,j,k)*Bprim(i,j,k,2)**2+g33(i,j,k)*Bprim(i,j,k,3)**2+ & + 2.0*(g12(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,2)+g13(i,j,k)*Bprim(i,j,k,1)*Bprim(i,j,k,3)+ & + g23(i,j,k)*Bprim(i,j,k,2)*Bprim(i,j,k,3)) + + dens(i,j,k) = sdet*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + rho(i,j,k) = GRHydro_rho_min + scon(i,j,k,:) = 0.d0 + vup(i,j,k,:) = 0.d0 + w_lorentz(i,j,k) = 1.d0 + + if(evolve_temper.ne.0) then + ! set hot atmosphere values + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + y_e_con(i,j,k) = y_e(i,j,k) * dens(i,j,k) + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr) + else + keytemp = 0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + endif + + ! w_lorentz=1, so the expression for tau reduces to: + + !!$ tau does need to take into account the existing B-field + !!$ with w_lorentz=1, we find tau = sqrtdet*(rho (1+eps+b^2/2)) - dens [Press drops out] + tau(i,j,k) = sdet * (rho(i,j,k)*(1.0+eps(i,j,k)+b2/2.0)) - dens(i,j,k) + + if(tau(i,j,k).le.sdet*b2*0.5d0)then + tau(i,j,k) = GRHydro_tau_min + sdet*b2*0.5d0 + endif + + end if + + end do + end do + end do + !$OMP END PARALLEL DO + + return +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_betax +#undef faulty_betay +#undef faulty_betaz +#undef faulty_vel +#undef faulty_Bvec + +end subroutine Conservative2PrimitiveAM + + + /*@@ + @routine Conservative2PrimitiveBoundariesAM + @date Sep 15, 2010 + @author Scott Noble, Joshua Faber, Bruno Mundim, The GRHydro Developers + @desc + This routine is used only if the reconstruction is performed on the conserved variables. + It computes the primitive variables on cell boundaries. + Since reconstruction on conservative had not proved to be very successful, + some of the improvements to the C2P routines (e.g. the check about + whether a failure happens in a point that will be restriced anyway) + are not implemented here yet. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + +subroutine Conservative2PrimitiveBoundsAM(CCTK_ARGUMENTS) + + use Con2PrimM_fortran_interfaces + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + integer :: i, j, k, itracer, nx, ny, nz + CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,& + uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin(1), epsmin(1) + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: b2minus, b2plus, local_gam(1), local_pgam,local_K,scminus,scplus + CCTK_INT :: epsnegative + CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp + character(len=100) warnline + + CCTK_REAL :: local_min_tracer + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress(1),xtemp(1),xye(1),xeps(1),xrho(1),one(1)=1.0d0 +! end EOS Omni vars + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define vel faulty_vel +#define Bvec faulty_Bvec + +! begin EOS Omni vars + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress(1)=0.0d0;xeps(1)=0.0d0;xtemp(1)=0.0d0;xye(1)=0.0d0 +! end EOS Omni vars + + ! this is a poly call + xrho(1)=GRHydro_rho_min + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,one,xtemp,xye,pmin,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,epsmin,xtemp,xye,pmin,epsmin,keyerr,anyerr) + + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + one,one,xtemp,xye,local_gam,keyerr,anyerr) + local_gam=local_gam+1.0 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + if (use_min_tracer .ne. 0) then + local_min_tracer = min_tracer + else + local_min_tracer = 0d0 + end if + + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + + !do not compute if in atmosphere or in an excised region + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & + GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle + + gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) + + epsnegative = 0 + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + +!!$ Tracers get no update for MHD! + if (evolve_tracer .ne. 0) then + do itracer=1,number_of_tracers + call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & + tracer(i,j,k,itracer), dens(i,j,k)) + enddo + + if (use_min_tracer .ne. 0) then + if (tracer(i,j,k,itracer) .le. local_min_tracer) then + tracer(i,j,k,itracer) = local_min_tracer + end if + end if + + endif + + if(evolve_Y_e.ne.0) then + Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) + endif + Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k) + Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k) + Bconsz_tmp = sqrt(avg_detl)*Bveczminus(i,j,k) + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densminus(i,j,k), & + sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, xye(1), xtemp(1), rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + if (epsnegative .ne. 0) then + !$OMP CRITICAL + call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!') + !$OMP END CRITICAL + + xrho=10.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + one,one,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress(1) + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,one,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress(1)/local_K)/log(xrho(1)) + scminus = local_K*densminus(i,j,k) + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densminus(i,j,k), & + sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus, & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + end if + + if (epsminus(i,j,k) .lt. 0.0d0) then + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + !$OMP CRITICAL + call CCTK_WARN(1,'Con2Prim: stopping the code.') + call CCTK_WARN(1, ' specific internal energy just went below 0! ') + write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,'(a20,3g16.7)') 'xyz location: ',& + x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,'(a20,3g16.7)') 'velocities: ',& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(GRHydro_c2p_warnlevel, "Specific internal energy negative") + !$OMP END CRITICAL + exit + endif + endif + + epsnegative = 0 + + Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k) + Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k) + Bconsz_tmp = sqrt(avg_detr)*Bveczplus(i,j,k) + call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, keytemp, GRHydro_eos_rf_prec, local_gam(1),densplus(i,j,k), & + sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k), & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, xye(1), xtemp(1), rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + if (epsnegative .ne. 0) then + !$OMP CRITICAL + call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!') + !$OMP END CRITICAL + + xrho=10.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + one,one,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress(1) + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,one,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress(1)/local_K)/log(xrho(1)) + scplus = local_K*densplus(i,j,k) + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densplus(i,j,k), & + sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end if + + if (epsplus(i,j,k) .lt. 0.0d0) then + if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then + !$OMP CRITICAL + call CCTK_WARN(1,'Con2Prim: stopping the code.') + call CCTK_WARN(1, ' specific internal energy just went below 0! ') + write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel + call CCTK_WARN(1,warnline) + write(warnline,'(a20,3g16.7)') 'xyz location: ',& + x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,'(a20,g16.7)') 'radius: ',r(i,j,k) + call CCTK_WARN(1,warnline) + write(warnline,'(a20,3g16.7)') 'velocities: ',& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k) + call CCTK_WARN(1,warnline) + call CCTK_WARN(GRHydro_c2p_warnlevel, "Specific internal energy negative") + write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',& + x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k) + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + endif + endif + end do + end do + end do + +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_betax +#undef faulty_betay +#undef faulty_betaz +#undef faulty_vel +#undef faulty_Bvec + +end subroutine Conservative2PrimitiveBoundsAM + +/*@@ +@routine Con2PrimPolytypeAM +@date Sep 16, 2010 +@author SCott Noble, Joshua Faber, Bruno Mundim, Ian Hawke +@desc +All routines below are identical to those above, just +specialised from polytropic type EOS. +@enddesc +@calls +@calledby +@history + +@endhistory + +@@*/ + +subroutine Conservative2PrimitivePolytypeAM(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k, itracer, nx, ny, nz + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det,b2 + + CCTK_INT :: epsnegative + CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp + + CCTK_REAL :: local_min_tracer, local_pgam,local_K, sc + ! character(len=400) :: warnline + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,xtemp,xye,xeps,xrho +! end EOS Omni vars + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + + call GRHydro_BvecfromAvec(CCTK_PASS_FTOF) + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define vel faulty_vel +#define Bvec faulty_Bvec + +! begin EOS Omni vars + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0 +! end EOS Omni vars + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + if (use_min_tracer .ne. 0) then + local_min_tracer = min_tracer + else + local_min_tracer = 0d0 + end if + +!!$ do k = GRHydro_stencil + 1, nz - GRHydro_stencil +!!$ do j = GRHydro_stencil + 1, ny - GRHydro_stencil +!!$ do i = GRHydro_stencil + 1, nx - GRHydro_stencil + + + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, & + !$OMP Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, & + !$OMP b2, xrho, xpress, local_K, local_pgam, sc) + do k = 1, nz + do j = 1, ny + do i = 1, nx + + !do not compute if in atmosphere or in an excised region + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & + GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle + + det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k)) + call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),& + g23(i,j,k),g33(i,j,k)) + +!!$ No MHD changes to tracers + if (evolve_tracer .ne. 0) then + do itracer=1,number_of_tracers + call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & + tracer(i,j,k,itracer), dens(i,j,k)) + enddo + + if (use_min_tracer .ne. 0) then + if (tracer(i,j,k,itracer) .le. local_min_tracer) then + tracer(i,j,k,itracer) = local_min_tracer + end if + end if + + endif + + if(evolve_Y_e.ne.0) then + Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k) + endif + + xrho=10.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + 1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress/local_K)/log(xrho) + sc = local_K*dens(i,j,k) + Bconsx_tmp = sqrt(det)*Bprim(i,j,k,1) + Bconsy_tmp = sqrt(det)*Bprim(i,j,k,2) + Bconsz_tmp = sqrt(det)*Bprim(i,j,k,3) + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & + scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rho(i,j,k),& + vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& + Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & + uxx,uxy,uxz,uyy,uyz,uzz,det, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + end do + end do + end do + + !$OMP END PARALLEL DO + + return + +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_betax +#undef faulty_betay +#undef faulty_betaz +#undef faulty_vel +#undef faulty_Bvec + +end subroutine Conservative2PrimitivePolytypeAM + + + /*@@ + @routine Cons2PrimBoundsPolytypeAM + @date Sep 16, 2010 + @author Scott Noble, Joshua Faber, Bruno Mundim, The GRHydro Developers + @desc + This routine is used only if the reconstruction is performed on the conserved variables. + It computes the primitive variables on cell boundaries. + Since reconstruction on conservative had not proved to be very successful, + some of the improvements to the C2P routines (e.g. the check about + whether a failure happens in a point that will be restriced anyway) are not implemented here yet. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Con2PrimBoundsPolytypeAM(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k, nx, ny, nz + CCTK_REAL :: uxxl, uxyl, uxzl, uyyl, uyzl, uzzl,& + uxxr, uxyr, uxzr, uyyr, uyzr, uzzr + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: b2minus, b2plus + CCTK_REAL :: Bconsx_tmp, Bconsy_tmp, Bconsz_tmp + CCTK_INT :: epsnegative + + CCTK_REAL :: local_pgam,local_K,scplus,scminus + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,xtemp,xye,xeps,xrho +! end EOS Omni vars + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define vel faulty_vel +#define Bvec faulty_Bvec + +! begin EOS Omni vars + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress=0.0d0;xtemp=0.0d0;xye=0.0d0;xeps=0.0d0 +! end EOS Omni vars + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 + do j = GRHydro_stencil, ny - GRHydro_stencil + 1 + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + + !do not compute if in atmosphere or in an excised region + if ((atmosphere_mask(i,j,k) .ne. 0) .or. & + GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle + + gxxl = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + + xrho=10.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + 1.d0,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress + + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress/local_K)/log(xrho) + scminus = local_K*densminus(i,j,k) + scplus = local_K*densplus(i,j,k) + Bconsx_tmp = sqrt(avg_detl)*Bvecxminus(i,j,k) + Bconsy_tmp = sqrt(avg_detl)*Bvecyminus(i,j,k) + Bconsz_tmp = sqrt(avg_detl)*Bveczminus(i,j,k) + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densminus(i,j,k), & + sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus,& + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),& + Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + + Bconsx_tmp = sqrt(avg_detr)*Bvecxplus(i,j,k) + Bconsy_tmp = sqrt(avg_detr)*Bvecyplus(i,j,k) + Bconsz_tmp = sqrt(avg_detr)*Bveczplus(i,j,k) + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), & + sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,& + Bconsx_tmp, Bconsy_tmp, Bconsz_tmp, rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,w_lorentzplus(i,j,k),& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, & + epsnegative,GRHydro_C2P_failed(i,j,k)) + end do + end do + end do + +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_betax +#undef faulty_betay +#undef faulty_betaz +#undef faulty_vel +#undef faulty_Bvec + +end subroutine Con2PrimBoundsPolytypeAM + +!!$ Con2Prim_ptTracer, Con2Prim_BoundsTracer, and Con2Prim_ptBoundsTracer need not be rewritten! + diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90 index 83372cb..c332aff 100644 --- a/src/GRHydro_Con2PrimM.F90 +++ b/src/GRHydro_Con2PrimM.F90 @@ -806,7 +806,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS) else local_min_tracer = 0d0 end if - + do k = GRHydro_stencil, nz - GRHydro_stencil + 1 do j = GRHydro_stencil, ny - GRHydro_stencil + 1 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 diff --git a/src/GRHydro_HLLE_AM.F90 b/src/GRHydro_HLLE_AM.F90 new file mode 100644 index 0000000..45e3afa --- /dev/null +++ b/src/GRHydro_HLLE_AM.F90 @@ -0,0 +1,943 @@ + /*@@ + @file GRHydro_HLLEPolyM.F90 + @date Aug 30, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke, Pedro Montero, Toni Font + @desc + The HLLE solver. Called from the wrapper function, so works in + all directions. + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#include "GRHydro_Macros.h" +#include "SpaceMask.h" + + /*@@ + @routine GRHydro_HLLE_AM + @date Aug 30, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke, Pedro Montero, Toni Font + @desc + The HLLE solver. Sufficiently simple that its just one big routine. + @enddesc + @calls + @calledby + @history + Altered from Cactus 3 routines originally written by Toni Font. + @endhistory + +@@*/ + +subroutine GRHydro_HLLE_AM(CCTK_ARGUMENTS) + USE GRHydro_EigenproblemM + USE GRHydro_Scalars + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k, m + CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff + CCTK_REAL, dimension(10) :: prim_p, prim_m + CCTK_REAL, dimension(5) :: lamminus,lamplus + CCTK_REAL :: charmin, charmax, charpm,chartop,avg_alp,avg_det, sdet + CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & + cs2_p, cs2_m, dpdeps_p, dpdeps_m + CCTK_REAL :: rhoenth_p, rhoenth_m, avg_betax, avg_betay, avg_betaz + CCTK_REAL :: vxtp,vytp,vztp,vxtm,vytm,vztm,ab0p,ab0m,b2p,b2m,bdotvp,bdotvm + CCTK_REAL :: wp,wm,v2p,v2m,bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,vA2m,vA2p + CCTK_REAL :: Bvecxlowp,Bvecxlowm,Bvecylowp,Bvecylowm,Bveczlowp,Bveczlowm + CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm + + CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm + CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc + + CCTK_INT :: type_bits, trivial, not_trivial + CCTK_REAL :: xtemp + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define vel faulty_vel +#define Bvec faulty_Bvec + + if (flux_direction == 1) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & + &"trivial") + else if (flux_direction == 2) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", & + &"trivial") + else if (flux_direction == 3) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & + &"trivial") + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction + + !$OMP PARALLEL DO PRIVATE(k,j,i,f1,lamminus,lamplus,cons_p,cons_m,fplus,fminus,qdiff,psidcf,psidcp,psidcm,prim_p,prim_m,& + !$OMP avg_betax,avg_betay,avg_betaz,avg_beta,avg_alp,& + !$OMP gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det,sdet,uxxh,uxyh,uxzh,uyyh,uyzh,uzzh,& + !$OMP vxtp,vxtm,vytp,vytm,vztp,vztm,& + !$OMP velxlowp,velxlowm,Bvecxlowp,Bvecxlowm,& + !$OMP velylowp,velylowm,Bvecylowp,Bvecylowm,& + !$OMP velzlowp,velzlowm,Bveczlowp,Bveczlowm,& + !$OMP bdotvp,bdotvm,b2p,b2m,v2p,v2m,wp,wm,& + !$OMP bxlowp,bxlowm,bylowp,bylowm,bzlowp,bzlowm,& + !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,usendh,psidcdiff,psidcfp,psidcfm,charmin,charmax,chartop,charpm,m,xtemp) + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + transport_constraints*(1-xoffset) + + f1 = 0.d0 + lamminus = 0.d0 + lamplus = 0.d0 + cons_p = 0.d0 + cons_m = 0.d0 + fplus = 0.d0 + fminus = 0.d0 + qdiff = 0.d0 + + if(clean_divergence.ne.0) then + psidcp = 0.d0 + psidcm = 0.d0 + endif + +!!$ Set the left (p for plus) and right (m_i for minus, i+1) states + + cons_p(1) = densplus(i,j,k) + cons_p(2) = sxplus(i,j,k) + cons_p(3) = syplus(i,j,k) + cons_p(4) = szplus(i,j,k) + cons_p(5) = tauplus(i,j,k) + cons_p(6) = Avecxplus(i,j,k) + cons_p(7) = Avecyplus(i,j,k) + cons_p(8) = Aveczplus(i,j,k) + + cons_m(1) = densminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(2) = sxminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(3) = syminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(4) = szminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(5) = tauminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(6) = Avecxminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(7) = Avecyminus(i+xoffset,j+yoffset,k+zoffset) + cons_m(8) = Aveczminus(i+xoffset,j+yoffset,k+zoffset) + + prim_p(1) = rhoplus(i,j,k) + prim_p(2) = velxplus(i,j,k) + prim_p(3) = velyplus(i,j,k) + prim_p(4) = velzplus(i,j,k) + prim_p(5) = epsplus(i,j,k) + prim_p(6) = pressplus(i,j,k) + prim_p(7) = w_lorentzplus(i,j,k) + prim_p(8) = Bvecxplus(i,j,k) + prim_p(9) = Bvecyplus(i,j,k) + prim_p(10) = Bveczplus(i,j,k) + + prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) + prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(7) = w_lorentzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(8) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(9) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(10)= Bveczminus(i+xoffset,j+yoffset,k+zoffset) + + if(clean_divergence.ne.0) then + psidcp = psidcplus(i,j,k) + psidcm = psidcminus(i+xoffset,j+yoffset,k+zoffset) + endif + +!!$ Calculate various metric terms here. +!!$ Note also need the average of the lapse at the +!!$ left and right points. +!!$ +!!$ In MHD, we need all three shift components regardless of the flux direction + + avg_betax = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + & + beta1(i,j,k)) + avg_betay = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + & + beta2(i,j,k)) + avg_betaz = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & + beta3(i,j,k)) + if (flux_direction == 1) then + avg_beta=avg_betax + else if (flux_direction == 2) then + avg_beta=avg_betay + else if (flux_direction == 3) then + avg_beta=avg_betaz + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) + + gxxh = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k)) + gxyh = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k)) + gxzh = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k)) + gyyh = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k)) + gyzh = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k)) + gzzh = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(i,j,k)) + + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) + sdet = sqrt(avg_det) + + call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & + avg_det,gxxh, gxyh, gxzh, & + gyyh, gyzh, gzzh) + + + vxtp = prim_p(2)-avg_betax/avg_alp + vytp = prim_p(3)-avg_betay/avg_alp + vztp = prim_p(4)-avg_betaz/avg_alp + vxtm = prim_m(2)-avg_betax/avg_alp + vytm = prim_m(3)-avg_betay/avg_alp + vztm = prim_m(4)-avg_betaz/avg_alp + + call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & + prim_p(2),prim_p(3),prim_p(4),prim_p(8),prim_p(9),prim_p(10), & + velxlowp,velylowp,velzlowp,Bvecxlowp,Bvecylowp,Bveczlowp, & + bdotvp,b2p,v2p,wp,bxlowp,bylowp,bzlowp) + call calc_vlow_blow(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, & + prim_m(2),prim_m(3),prim_m(4),prim_m(8),prim_m(9),prim_m(10), & + velxlowm,velylowm,velzlowm,Bvecxlowm,Bvecylowm,Bveczlowm, & + bdotvm,b2m,v2m,wm,bxlowm,bylowm,bzlowm) + + rhoenth_p = prim_p(1)*(1.d0+prim_p(5))+prim_p(6) + rhoenth_m = prim_m(1)*(1.d0+prim_m(5))+prim_m(6) + + ab0p = wp*bdotvp + ab0m = wm*bdotvm + + vA2p = b2p/(rhoenth_p+b2p) + vA2m = b2m/(rhoenth_m+b2m) + +!!$ p^* = p+b^2/2 in Anton et al. + pressstarp = prim_p(6)+0.5d0*b2p + pressstarm = prim_m(6)+0.5d0*b2m + + +!!$ If the Riemann problem is trivial, just calculate the fluxes from the +!!$ left state and skip to the next cell + + if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then + +!!$ we now pass in the B-field as conserved and flux, the vtilde's instead of v's, +!!$ pressstar instead of P, b_i, alp b^0, w, metric determinant, +!!$ alp, and beta in the flux dir + + if (flux_direction == 1) then + call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& + cons_p(6),cons_p(7),cons_p(8),& + f1(1),f1(2),f1(3),f1(4),f1(5),f1(6),f1(7),f1(8),& + vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + if(clean_divergence.ne.0) then + f1(6)=f1(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + psidcf = cons_p(6)/sdet-psidcp*avg_betax/avg_alp + endif + + else if (flux_direction == 2) then + call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& + cons_p(7),cons_p(8),cons_p(6),& + f1(1),f1(3),f1(4),f1(2),f1(5),f1(7),f1(8),f1(6),& + vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + if(clean_divergence.ne.0) then + f1(6)=f1(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + psidcf = cons_p(7)/sdet-psidcp*avg_betay/avg_alp + endif + + else if (flux_direction == 3) then + call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& + cons_p(8),cons_p(6),cons_p(7),& + f1(1),f1(4),f1(2),f1(3),f1(5),f1(8),f1(6),f1(7), & + vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + if(clean_divergence.ne.0) then + f1(6)=f1(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp + f1(7)=f1(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp + f1(8)=f1(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + psidcf = cons_p(8)/sdet-psidcp*avg_betaz/avg_alp + endif + + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + else !!! The end of this branch is right at the bottom of the routine + + if (flux_direction == 1) then + usendh = uxxh + else if (flux_direction == 2) then + usendh = uyyh + else if (flux_direction == 3) then + usendh = uzzh + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + +!!$ Calculate the jumps in the conserved variables + + qdiff(1) = cons_m(1) - cons_p(1) + qdiff(2) = cons_m(2) - cons_p(2) + qdiff(3) = cons_m(3) - cons_p(3) + qdiff(4) = cons_m(4) - cons_p(4) + qdiff(5) = cons_m(5) - cons_p(5) + qdiff(6) = cons_m(6) - cons_p(6) + qdiff(7) = cons_m(7) - cons_p(7) + qdiff(8) = cons_m(8) - cons_p(8) + + if (clean_divergence.ne.0) then + psidcdiff = psidcm - psidcp + endif + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + if(evolve_temper.ne.1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(8),prim_m(9),prim_m(10),& + lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(8),prim_p(9),prim_p(10),& + lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + prim_m(8),prim_m(9),prim_m(10),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + prim_p(8),prim_p(9),prim_p(10),& + xtemp,y_e_plus(i,j,k),& + lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + endif + + call num_x_fluxM(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& + cons_p(6),cons_p(7),cons_p(8),& + fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),fplus(6),fplus(7),fplus(8),& + vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + call num_x_fluxM(cons_m(1),cons_m(2),cons_m(3),cons_m(4),cons_m(5),& + cons_m(6),cons_m(7),cons_m(8),& + fminus(1),fminus(2),fminus(3),fminus(4),fminus(5),& + fminus(6),fminus(7),fminus(8),& + vxtm,vytm,vztm,pressstarm,bxlowm,bylowm,bzlowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + + if(clean_divergence.ne.0) then + fminus(6)=fminus(6) + 1.0d0*sdet*uxxh*psidcm - cons_m(6)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0d0*sdet*uxyh*psidcm - cons_m(6)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0d0*sdet*uxzh*psidcm - cons_m(6)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0d0*sdet*uxxh*psidcp - cons_p(6)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0d0*sdet*uxyh*psidcp - cons_p(6)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0d0*sdet*uxzh*psidcp - cons_p(6)*avg_betaz/avg_alp + psidcfp = cons_p(6)/sdet-avg_betax*psidcp/avg_alp + psidcfm = cons_m(6)/sdet-avg_betax*psidcm/avg_alp + endif + + else if (flux_direction == 2) then + if(evolve_temper.ne.1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(9),prim_m(10),prim_m(8),& + lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(9),prim_p(10),prim_p(8),& + lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + prim_m(9),prim_m(10),prim_m(8),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + prim_p(9),prim_p(10),prim_p(8),& + xtemp,y_e_plus(i,j,k),& + lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + endif + + call num_x_fluxM(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& + cons_p(7),cons_p(8),cons_p(6),& + fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),fplus(7),fplus(8),fplus(6),& + vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + call num_x_fluxM(cons_m(1),cons_m(3),cons_m(4),cons_m(2),cons_m(5),& + cons_m(7),cons_m(8),cons_m(6),& + fminus(1),fminus(3),fminus(4),fminus(2),fminus(5),& + fminus(7),fminus(8),fminus(6),& + vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + + if(clean_divergence.ne.0) then + fminus(6)=fminus(6) + 1.0d0*sdet*uxyh*psidcm - cons_m(7)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0d0*sdet*uyyh*psidcm - cons_m(7)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0d0*sdet*uyzh*psidcm - cons_m(7)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0d0*sdet*uxyh*psidcp - cons_p(7)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0d0*sdet*uyyh*psidcp - cons_p(7)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0d0*sdet*uyzh*psidcp - cons_p(7)*avg_betaz/avg_alp + psidcfp = cons_p(7)/sdet-avg_betay*psidcp/avg_alp + psidcfm = cons_m(7)/sdet-avg_betay*psidcm/avg_alp + endif + + else if (flux_direction == 3) then + if(evolve_temper.ne.1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(10),prim_m(8),prim_m(9),& + lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(10),prim_p(8),prim_p(9),& + lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + else + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + prim_m(10),prim_m(8),prim_m(9),& + xtemp,y_e_minus(i+xoffset,j+yoffset,k+zoffset),& + lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + xtemp = temperature(i,j,k) + call eigenvaluesM_hot(GRHydro_eos_handle,i,j,k,& + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + prim_p(10),prim_p(8),prim_p(9),& + xtemp,y_e_plus(i,j,k),& + lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + endif + + call num_x_fluxM(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& + cons_p(8),cons_p(6),cons_p(7),& + fplus(1),fplus(4),fplus(2),fplus(3),fplus(5),fplus(8),fplus(6),fplus(7), & + vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, & + avg_det,avg_alp,avg_beta) + call num_x_fluxM(cons_m(1),cons_m(4),cons_m(2),cons_m(3),cons_m(5),& + cons_m(8),cons_m(6),cons_m(7),& + fminus(1),fminus(4),fminus(2),fminus(3),fminus(5), & + fminus(8),fminus(6),fminus(7), & + vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, & + avg_det,avg_alp,avg_beta) + + if(clean_divergence.ne.0) then + fminus(6)=fminus(6) + 1.0d0*sdet*uxzh*psidcm - cons_m(8)*avg_betax/avg_alp + fminus(7)=fminus(7) + 1.0d0*sdet*uyzh*psidcm - cons_m(8)*avg_betay/avg_alp + fminus(8)=fminus(8) + 1.0d0*sdet*uzzh*psidcm - cons_m(8)*avg_betaz/avg_alp + fplus(6)=fplus(6) + 1.0d0*sdet*uxzh*psidcp - cons_p(8)*avg_betax/avg_alp + fplus(7)=fplus(7) + 1.0d0*sdet*uyzh*psidcp - cons_p(8)*avg_betay/avg_alp + fplus(8)=fplus(8) + 1.0d0*sdet*uzzh*psidcp - cons_p(8)*avg_betaz/avg_alp + psidcfp = cons_p(8)/sdet-avg_betaz*psidcp/avg_alp + psidcfm = cons_m(8)/sdet-avg_betaz*psidcm/avg_alp + endif + + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + +!!$ Find minimum and maximum wavespeeds + + charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + chartop = max(-charmin,charmax) + + charpm = charmax - charmin + +!!$ Calculate flux by standard formula + + do m = 1,8 + + qdiff(m) = cons_m(m) - cons_p(m) + + if (HLLE) then + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + else if (LLF) then + f1(m) = 0.5d0 * (fplus(m) + fminus(m) - chartop * qdiff(m)) + end if + + end do + + if(clean_divergence.ne.0) then + + psidcdiff = psidcm - psidcp + select case(whichpsidcspeed) + case(0) + if (HLLE) then + psidcf = (charmax * psidcfp - charmin * psidcfm + & + charmax * charmin * psidcdiff) / charpm + else if (LLF) then + psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff) + end if + case(1) + !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? + !!$ psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & + !!$ 1.d0 * (-1.d0) * psidcdiff) + + !!$ The fastest speed for psidc comes from the condition + !!$ that the normal vector to the characteristic hypersurface + !!$ be spacelike (Eq. 60 of Anton et al.) + + charmax_dc = sqrt(usendh) - avg_beta/avg_alp + charmin_dc = -1.d0*sqrt(usendh) - avg_beta/avg_alp + charpm_dc = charmax_dc - charmin_dc + + psidcf = (charmax_dc * psidcfp - charmin_dc * psidcfm + & + charmax_dc * charmin_dc * psidcdiff) / charpm_dc + + if(decouple_normal_Bfield .ne. 0) then ! same expression for HLLE and LLF + !!$ B^i field decouples from the others and has same propagation + !!$ speed as divergence -null direction, + !!$ \pm sqrt(g^{xx}} - beta^x/alpha + f1(5+flux_direction) = (charmax_dc * fplus(5+flux_direction) & + - charmin_dc * fminus(5+flux_direction) + & + charmax_dc * charmin_dc * qdiff(5+flux_direction)) / charpm_dc + end if + case(2) + charmax = setcharmax + charmin = setcharmin + if (HLLE) then + psidcf = (charmax * psidcfp - charmin * psidcfm + & + charmax * charmin * psidcdiff) / charpm + else if (LLF) then + chartop = max(-charmin,charmax) + psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff) + end if + end select + endif + + + end if !!! The end of the SpaceMask check for a trivial RP. + + densflux(i, j, k) = f1(1) + sxflux(i, j, k) = f1(2) + syflux(i, j, k) = f1(3) + szflux(i, j, k) = f1(4) + tauflux(i, j, k) = f1(5) + + if ( evolve_Lorenz_gge.gt.0 ) then + !! FIX: These aren't zero + Avecxflux(i,j,k) = 0.d0 + Avecyflux(i,j,k) = 0.d0 + Aveczflux(i,j,k) = 0.d0 + Aphiflux(i,j,k) = 0.d0 + else + Avecxflux(i,j,k) = 0.d0 + Avecyflux(i,j,k) = 0.d0 + Aveczflux(i,j,k) = 0.d0 + end if + + if(clean_divergence.ne.0) then + psidcflux(i,j,k) = psidcf + endif + + if(evolve_Y_e.ne.0) then + if (densflux(i, j, k) > 0.d0) then + Y_e_con_flux(i, j, k) = & + Y_e_plus(i, j, k) * & + densflux(i, j, k) + else + Y_e_con_flux(i, j, k) = & + Y_e_minus(i + xoffset, j + yoffset, k + zoffset) * & + densflux(i, j, k) + endif + endif + + end do + end do + end do + !$OMP END PARALLEL DO +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_betax +#undef faulty_betay +#undef faulty_betaz +#undef faulty_vel +#undef faulty_Bvec + +end subroutine GRHydro_HLLE_AM + + /*@@ + @routine GRHydro_HLLE_TracerAM + @date Aug 30, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + HLLE just for the tracer. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_HLLE_TracerAM(CCTK_ARGUMENTS) + + USE GRHydro_EigenproblemM + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k, m + CCTK_REAL, dimension(number_of_tracers) :: cons_p,cons_m,fplus,fminus,f1 + CCTK_REAL, dimension(5) :: lamminus,lamplus + CCTK_REAL, dimension(number_of_tracers) :: qdiff + CCTK_REAL, dimension(7) :: prim_p, prim_m + CCTK_REAL, dimension(3) :: mag_p, mag_m + CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det + CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & + cs2_p, cs2_m, dpdeps_p, dpdeps_m + CCTK_REAL :: b2p,b2m,vA2m,vA2p + + CCTK_INT :: type_bits, trivial, not_trivial + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + beta1 => betaa + beta2 => betab + beta3 => betac + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + beta1 => betax + beta2 => betay + beta3 => betaz + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define vel faulty_vel +#define Bvec faulty_Bvec + + if (flux_direction == 1) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemX") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemX", & + &"trivial") + else if (flux_direction == 2) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemY") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemY", & + &"trivial") + else if (flux_direction == 3) then + call SpaceMask_GetTypeBits(type_bits, "Hydro_RiemannProblemZ") + call SpaceMask_GetStateBits(trivial, "Hydro_RiemannProblemZ", & + &"trivial") + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + + f1 = 0.d0 + lamminus = 0.d0 + lamplus = 0.d0 + cons_p = 0.d0 + cons_m = 0.d0 + mag_p = 0.d0 + mag_m = 0.d0 + fplus = 0.d0 + fminus = 0.d0 + qdiff = 0.d0 + +!!$ Set the left (p for plus) and right (m_i for minus, i+1) states + + cons_p(:) = cons_tracerplus(i,j,k,:) + cons_m(:) = cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + + mag_p(1) = Bvecxplus(i,j,k) + mag_p(2) = Bvecyplus(i,j,k) + mag_p(3) = Bveczplus(i,j,k) + + mag_m(1) = Bvecxminus(i+xoffset,j+yoffset,k+zoffset) + mag_m(2) = Bvecyminus(i+xoffset,j+yoffset,k+zoffset) + mag_m(3) = Bveczminus(i+xoffset,j+yoffset,k+zoffset) + + prim_p(1) = rhoplus(i,j,k) + prim_p(2) = velxplus(i,j,k) + prim_p(3) = velyplus(i,j,k) + prim_p(4) = velzplus(i,j,k) + prim_p(5) = epsplus(i,j,k) + prim_p(6) = pressplus(i,j,k) + prim_p(7) = w_lorentzplus(i,j,k) + + prim_m(1) = rhominus(i+xoffset,j+yoffset,k+zoffset) + prim_m(2) = velxminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(3) = velyminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(4) = velzminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(5) = epsminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(6) = pressminus(i+xoffset,j+yoffset,k+zoffset) + prim_m(7) = w_lorentzminus(i+xoffset,j+yoffset,k+zoffset) + +!!$ Calculate various metric terms here. +!!$ Note also need the average of the lapse at the +!!$ left and right points. + + if (flux_direction == 1) then + avg_beta = 0.5d0 * (beta1(i+xoffset,j+yoffset,k+zoffset) + & + beta1(i,j,k)) + else if (flux_direction == 2) then + avg_beta = 0.5d0 * (beta2(i+xoffset,j+yoffset,k+zoffset) + & + beta2(i,j,k)) + else if (flux_direction == 3) then + avg_beta = 0.5d0 * (beta3(i+xoffset,j+yoffset,k+zoffset) + & + beta3(i,j,k)) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + + avg_alp = 0.5 * (alp(i,j,k) + alp(i+xoffset,j+yoffset,k+zoffset)) + + gxxh = 0.5d0 * (g11(i+xoffset,j+yoffset,k+zoffset) + g11(i,j,k)) + gxyh = 0.5d0 * (g12(i+xoffset,j+yoffset,k+zoffset) + g12(i,j,k)) + gxzh = 0.5d0 * (g13(i+xoffset,j+yoffset,k+zoffset) + g13(i,j,k)) + gyyh = 0.5d0 * (g22(i+xoffset,j+yoffset,k+zoffset) + g22(i,j,k)) + gyzh = 0.5d0 * (g23(i+xoffset,j+yoffset,k+zoffset) + g23(i,j,k)) + gzzh = 0.5d0 * (g33(i+xoffset,j+yoffset,k+zoffset) + g33(i,j,k)) + + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) + + call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & + avg_det,gxxh, gxyh, gxzh, & + gyyh, gyzh, gzzh) + + if (flux_direction == 1) then + usendh = uxxh + else if (flux_direction == 2) then + usendh = uyyh + else if (flux_direction == 3) then + usendh = uzzh + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + +!!$ b^2 = (1-v^2)B^2+(B dot v)^2 + b2p=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_p(1),mag_p(2),mag_p(3))/prim_p(7)**2 + & + (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_p(2),prim_p(3),prim_p(4),mag_p(1),mag_p(2),mag_p(3)))**2 + b2m=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3))/prim_m(7)**2 + & + (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3)))**2 + + vA2p = b2p/(prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)+b2p) + vA2m = b2m/(prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)+b2m) + +!!$ Calculate the jumps in the conserved variables + + qdiff = cons_m - cons_p + +!!$ Eigenvalues and fluxes either side of the cell interface + + if (flux_direction == 1) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(2),prim_m(3),prim_m(4),prim_m(5),prim_m(7), & + mag_m(1),mag_m(2),mag_m(3),& + lamminus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(2),prim_p(3),prim_p(4),prim_p(5),prim_p(7), & + mag_p(1),mag_p(2),mag_p(3),& + lamplus,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,& + usendh,avg_alp,avg_beta) + fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * & + cons_tracerplus(i,j,k,:) + fminus(:) = (velxminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & + cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + else if (flux_direction == 2) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(3),prim_m(4),prim_m(2),prim_m(5),prim_m(7), & + mag_m(2),mag_m(3),mag_m(1),& + lamminus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle, & + prim_p(1),prim_p(3),prim_p(4),prim_p(2),prim_p(5),prim_p(7), & + mag_p(2),mag_p(3),mag_p(1),& + lamplus,gyyh,gyzh,gxyh,gzzh,gxzh,gxxh,& + usendh,avg_alp,avg_beta) + fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * & + cons_tracerplus(i,j,k,:) + fminus(:) = (velyminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & + cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + else if (flux_direction == 3) then + call eigenvaluesM(GRHydro_eos_handle,& + prim_m(1),prim_m(4),prim_m(2),prim_m(3),prim_m(5),prim_m(7), & + mag_m(3),mag_m(1),mag_m(2),& + lamminus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + call eigenvaluesM(GRHydro_eos_handle,& + prim_p(1),prim_p(4),prim_p(2),prim_p(3),prim_p(5),prim_p(7), & + mag_p(3),mag_p(1),mag_p(2),& + lamplus,gzzh,gxzh,gyzh,gxxh,gxyh,gyyh,& + usendh,avg_alp,avg_beta) + fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * & + cons_tracerplus(i,j,k,:) + fminus(:) = (velzminus(i+xoffset,j+yoffset,k+zoffset) - avg_beta / avg_alp) * & + cons_tracerminus(i+xoffset,j+yoffset,k+zoffset,:) + else + call CCTK_WARN(0, "Flux direction not x,y,z") + end if + +!!$ Find minimum and maximum wavespeeds + + charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & + lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& + lamminus(4),lamminus(5)) + + charpm = charmax - charmin + +!!$ Calculate flux by standard formula + + do m = 1,number_of_tracers + + qdiff(m) = cons_m(m) - cons_p(m) + + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + + end do + + cons_tracerflux(i, j, k,:) = f1(:) +!!$ +!!$ if ( ((flux_direction.eq.3).and.(i.eq.4).and.(j.eq.4)).or.& +!!$ ((flux_direction.eq.2).and.(i.eq.4).and.(k.eq.4)).or.& +!!$ ((flux_direction.eq.1).and.(j.eq.4).and.(k.eq.4))& +!!$ ) then +!!$ write(*,*) flux_direction, i, j, k, f1(1), cons_m(1), cons_p(1) +!!$ end if + + end do + end do +end do +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_betax +#undef faulty_betay +#undef faulty_betaz +#undef faulty_vel +#undef faulty_Bvec + + +end subroutine GRHydro_HLLE_TracerAM + diff --git a/src/GRHydro_InterfacesAM.h b/src/GRHydro_InterfacesAM.h new file mode 100644 index 0000000..c2565d8 --- /dev/null +++ b/src/GRHydro_InterfacesAM.h @@ -0,0 +1,95 @@ +module Con2PrimAM_fortran_interfaces + implicit none + + interface + + subroutine GRHydro_Con2PrimM_pt( handle, keytemp, prec, & + local_gam, dens, & + sx, sy, sz, & + tau, & + Bconsx, Bconsy, Bconsz, & + xye, xtemp, & + rho, & + velx, vely, velz,& + epsilon, pressure,& + Bx, By, Bz, & + bsq,& + w_lorentz, & + gxx, gxy, gxz, & + gyy, gyz, gzz, & + uxx, uxy, uxz,& + uyy, uyz, uzz,& + det,& + epsnegative, & + retval) + + implicit none + CCTK_INT handle + CCTK_INT keytemp + CCTK_REAL prec + CCTK_REAL local_gam + CCTK_REAL dens + CCTK_REAL sx, sy, sz + CCTK_REAL tau + CCTK_REAL Bconsx, Bconsy, Bconsz + CCTK_REAL xye + CCTK_REAL xtemp + CCTK_REAL rho + CCTK_REAL velx, vely, velz + CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq + CCTK_REAL w_lorentz + CCTK_REAL gxx, gxy, gxz + CCTK_REAL gyy, gyz, gzz + CCTK_REAL uxx, uxy, uxz + CCTK_REAL uyy, uyz, uzz + CCTK_REAL det + CCTK_INT epsnegative + CCTK_REAL retval + end subroutine GRHydro_Con2PrimM_pt + + subroutine GRHydro_Con2PrimM_Polytype_pt( handle, local_gam,& + dens, & + sx, sy, sz, & + sc, & + Bconsx, Bconsy, Bconsz, & + rho, & + velx, vely, velz,& + epsilon, pressure,& + Bx, By, Bz, & + bsq,& + w_lorentz, & + gxx, gxy, gxz, & + gyy, gyz, gzz, & + uxx, uxy, uxz,& + uyy, uyz, uzz,& + det,& + epsnegative, & + retval) + + implicit none + CCTK_INT handle + CCTK_REAL local_gam + CCTK_REAL dens + CCTK_REAL sx, sy, sz + CCTK_REAL sc + CCTK_REAL Bconsx, Bconsy, Bconsz + CCTK_REAL rho + CCTK_REAL velx, vely, velz + CCTK_REAL epsilon, pressure + CCTK_REAL Bx, By, Bz + CCTK_REAL bsq + CCTK_REAL w_lorentz + CCTK_REAL gxx, gxy, gxz + CCTK_REAL gyy, gyz, gzz + CCTK_REAL uxx, uxy, uxz + CCTK_REAL uyy, uyz, uzz + CCTK_REAL det + CCTK_INT epsnegative + CCTK_REAL retval + end subroutine GRHydro_Con2PrimM_Polytype_pt + + end interface + +end module Con2PrimAM_fortran_interfaces diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90 index c738f79..1ab2a40 100644 --- a/src/GRHydro_ParamCheck.F90 +++ b/src/GRHydro_ParamCheck.F90 @@ -88,12 +88,20 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) call CCTK_PARAMWARN("GRHydro::bound = 'static' is no longer supported, use 'none' instead"); end if - if (CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then + if (CCTK_EQUALS(Bvec_evolution_method,"GRHydro").or.CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec")) then evolve_MHD = 1 else evolve_MHD = 0 endif + if ( CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") ) then + if ( CCTK_EQUALS(Avec_gauge,"Lorenz") ) then + evolve_Lorenz_gge = 1 + else + evolve_Lorenz_gge = 0 + endif + endif + if (CCTK_EQUALS(Y_e_evolution_method,"GRHydro")) then evolve_Y_e = 1 else diff --git a/src/GRHydro_Prim2ConAM.F90 b/src/GRHydro_Prim2ConAM.F90 new file mode 100644 index 0000000..35b2619 --- /dev/null +++ b/src/GRHydro_Prim2ConAM.F90 @@ -0,0 +1,762 @@ + /*@@ + @file primitive2conservative + @date Aug 31, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke + @desc + Primitive to conservative routine + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Avecx(i,j,k) Avec(i,j,k,1) +#define Avecy(i,j,k) Avec(i,j,k,2) +#define Avecz(i,j,k) Avec(i,j,k,3) + +#define DOT(x1,y1,z1,x2,y2,z2) ( DOTP(gxx,gxy,gxz,gyy,gyz,gzz,x1,y1,z1,x2,y2,z2) ) +#define DOT2(x1,y1,z1) ( DOTP2(gxx,gxy,gxz,gyy,gyz,gzz,x1,y1,z1) ) +#define DOTPT(x1,y1,z1,x2,y2,z2) ( DOTP(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt,x1,y1,z1,x2,y2,z2) ) +#define DOTPT2(x1,y1,z1) ( DOTP2(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt,x1,y1,z1) ) + + /*@@ + @routine primitive2conservativeAM + @date Aug 31 + @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke + @desc + Converts primitive to conserved variables for the boundary extended data. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine primitive2conservativeAM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + !CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,& + ! g11r,g12r,g13r,g22r,g23r,g33r + CCTK_REAL :: xtemp(1) + character(len=256) NaN_WarnLine + !TARGET gxx, gxy, gxz, gyy, gyz, gzz + + !CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + + !g11 => gxx + !g12 => gxy + !g13 => gxz + !g22 => gyy + !g23 => gyz + !g33 => gzz + + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + + if(evolve_temper.ne.1) then + + !$OMP PARALLEL DO PRIVATE(k,j,i,avg_detl,avg_detr,& + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,& + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) + + gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + + call prim2conAM(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& + rhominus(i,j,k), & + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & + Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k)) + + call prim2conAM(GRHydro_eos_handle, gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& + rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k), & + w_lorentzplus(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + + else + + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) + + gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + ! We will instead use an average temperature as an initial + ! guess. + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i-xoffset,j-yoffset,k-zoffset)) + + !$OMP CRITICAL + if (y_e_minus(i,j,k) .le. 0.0d0 .or. y_e_plus(i,j,k) .le. 0.0d0) then + write(NaN_WarnLine,'(a100,7g15.6)') '(y_e_minus,y_e_plus,x,y,z,rho)', y_e(i,j,k), y_e_minus(i,j,k), y_e_plus(i,j,k), x(i,j,k),y(i,j,k),z(i,j,k),rho(i,j,k) + call CCTK_WARN(1, NaN_WarnLine) + endif + !$OMP END CRITICAL + call prim2conAM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& + rhominus(i,j,k), & + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & + Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k), xtemp, y_e_minus(i,j,k)) + ! we do not have plus/minus vars for temperature since + ! eps is reconstructed. Hence, we do not update the + ! variable 'temperature' in prim2con at the interfaces + ! We will instead use an average temperature as an initial + ! guess. + xtemp(1) = 0.5d0*(temperature(i,j,k) + & + temperature(i+xoffset,j+yoffset,k+zoffset)) + + call prim2conAM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& + rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k), & + w_lorentzplus(i,j,k), xtemp, y_e_plus(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + endif + + +end subroutine primitive2conservativeAM + + /*@@ + @routine prim2conAM + @date Aug 31, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke + @desc + Converts from primitive to conservative at a single point + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine prim2conAM_hot(handle, GRHydro_reflevel, ii, jj, kk, x, y, z, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w, temp, ye) + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & + drho, dvelx, dvely, dvelz,& + deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz + CCTK_REAL :: temp(1),ye(1), yein(1), x, y, z + CCTK_INT :: handle, GRHydro_reflevel, ii, jj, kk + CCTK_REAL :: w + CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz + character(len=256) NaN_WarnLine + character(len=512) warnline + +! begin EOS Omni vars + integer :: n, keytemp, anyerr, keyerr(1) + ! real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) + real*8 :: temp0(1) + n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 + !xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 +! end EOS Omni vars + + temp0 = temp + w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1))) + +!!$ BEGIN: Check for NaN value + if (w .ne. w) then + !$OMP CRITICAL + write(NaN_WarnLine,'(a100,6g15.6)') 'NaN produced in sqrt(): (gxx,gxy,gxz,gyy,gyz,gzz)', gxx, gxy, gxz, gyy, gyz, gzz + call CCTK_WARN(1, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz + call CCTK_WARN(1, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (x,y,z)', x, y, z + call CCTK_WARN(1, NaN_WarnLine) + !write(NaN_WarnLine,'(a100,g15.6)') 'NaN produced in sqrt(): w', w + !call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) + !$OMP END CRITICAL + endif +!!$ END: Check for NaN value + + ye = max(min(ye,GRHydro_Y_e_max),GRHydro_Y_e_min) + + yein = ye + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + ! error handling + if(anyerr.ne.0) then + if(GRHydro_reflevel.lt.GRHydro_c2p_warn_from_reflevel) then + ! in this case (coarse grid error that is hopefully restricted + ! away), we use the average temperature between cells and call + ! the EOS with keytemp=1 + keytemp=1 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye,yein + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + endif + else + ! This is a way of recovering even on finer refinement levels: + ! Use the average temperature at the interface instead of the + ! reconstructed specific internal energy. + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp0,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + keytemp=1 + temp = temp0 + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,temp,ye,dpress,keyerr,anyerr) + keytemp=0 + if(anyerr.ne.0) then + !$OMP CRITICAL + call CCTK_WARN(1,"EOS error in prim2con_hot") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(1,warnline) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + !$OMP END CRITICAL + endif + endif + endif + + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz + vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz + vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz + + Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) + ab0=w*Bdotv + b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 + blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & + w*Bdotv*vlowx + blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & + w*Bdotv*vlowy + blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & + w*Bdotv*vlowz + + ddens = sqrt(det) * drho * w + dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ab0*blowx) + dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + ab0*blowy) + dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + ab0*blowz) + dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + + !!$OMP CRITICAL + !write(NaN_WarnLine,'(a100,3g15.6)') '(dens out, sqrt(det), w:)', ddens,sqrt(det),w + !call CCTK_WARN(1, NaN_WarnLine) + !!$OMP END CRITICAL +end subroutine prim2conAM_hot + + +subroutine prim2conAM(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL, dimension(1) :: ddens, dsx, dsy, dsz, dtau, & + drho, dvelx, dvely, dvelz,& + deps, dpress, dBvcx, dBvcy, dBvcz, vlowx, vlowy, vlowz + CCTK_REAL :: w + CCTK_REAL, dimension(1) :: Bdotv,ab0,b2,blowx,blowy,blowz + CCTK_INT :: handle + character(len=256) NaN_WarnLine + +! begin EOS Omni vars + integer :: n, keytemp, anyerr, keyerr(1) + real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) + n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 + xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 +! end EOS Omni vars + + w = 1.d0 / sqrt(1.d0 - DOT2(dvelx(1),dvely(1),dvelz(1))) + +!!$ BEGIN: Check for NaN value + if (w .ne. w) then + !$OMP CRITICAL + write(NaN_WarnLine,'(a100,6g15.6)') 'NaN produced in sqrt(): (gxx,gxy,gxz,gyy,gyz,gzz)', gxx, gxy, gxz, gyy, gyz, gzz + call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) + write(NaN_WarnLine,'(a100,3g15.6)') 'NaN produced in sqrt(): (dvelx,dvely,dvelz)', dvelx, dvely, dvelz + call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) + !$OMP END CRITICAL + endif +!!$ END: Check for NaN value + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,deps,xtemp,xye,dpress,keyerr,anyerr) + + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz + vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz + vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz + + Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) + ab0=w*Bdotv + b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 + blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & + w*Bdotv*vlowx + blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & + w*Bdotv*vlowy + blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & + w*Bdotv*vlowz + + ddens = sqrt(det) * drho * w + dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ab0*blowx) + dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + ab0*blowy) + dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + ab0*blowz) + dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + +end subroutine prim2conAM + + + /*@@ + @routine Primitive2ConservativeCellsAM + @date Aug 31, 2010 + @author + @desc + Wrapper function that converts primitive to conservative at the + cell centres. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine Primitive2ConservativeCellsAM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_REAL :: xtemp(1) + CCTK_INT :: i, j, k + CCTK_REAL :: det + CCTK_REAL :: maxtau0 + + maxtau0 = -1.0d60 + + if(evolve_temper.ne.1) then + + !$OMP PARALLEL DO PRIVATE(k,j,i,det), REDUCTION(MAX:maxtau0) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call prim2conAM(GRHydro_eos_handle,gxx(i,j,k),& + gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),& + rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & + Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k)) + + maxtau0 = max(maxtau0,tau(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + + ! TODO: to actually reduce GRHydro_tau_min over all Carpet components + ! we need to modify Carpet looping to allow this function to be called + ! on all AMR levels before calling any other function. The best would be + ! to set a special bin where functions would be called on all levels first + ! instead of calling all functions per level. The workaround for this problem + ! is to set GRHydro_tau_min to a user specified value as it was set in + ! GRHydro_Minima.F90. Once this issue is solved, then uncomment the line + ! below and create two other routines to be run in global mode so that + ! GRHydro_tau_min can be properly initialized and reduced. + + !GRHydro_tau_min = GRHydro_tau_min * maxtau0 + + else + !$OMP PARALLEL DO PRIVATE(k,j,i,det) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call prim2conAM_hot(GRHydro_eos_handle,GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),& + rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),Bvecx(i,j,k), & + Bvecy(i,j,k), Bvecz(i,j,k), w_lorentz(i,j,k), temperature(i,j,k), Y_e(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + endif + + if(evolve_Y_e.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k) + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif + + + +end subroutine Primitive2ConservativeCellsAM + + + /*@@ + @routine Prim2ConservativePolytypeAM + @date Aug 31, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + Same as first routine, only for polytropes. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + +subroutine Prim2ConservativePolytypeAM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer :: i, j, k + CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + + ! constraint transport needs to be able to average fluxes in the directions + ! other that flux_direction, which in turn need the primitives on interfaces + !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& + !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) + do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1 + transport_constraints*(1-zoffset) + do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1 + transport_constraints*(1-yoffset) + do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1 + transport_constraints*(1-xoffset) + + gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) + gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) + gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) + gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) + gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) + gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) + gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) + gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) + gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) + gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) + gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) + gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + + call prim2conpolytypeAM(GRHydro_eos_handle, gxxl,gxyl,gxzl,& + gyyl,gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& + rhominus(i,j,k),& + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & + Bvecyminus(i,j,k),Bveczminus(i,j,k),w_lorentzminus(i, j, k)) + + call prim2conpolytypeAM(GRHydro_eos_handle, gxxr,gxyr,gxzr,& + gyyr,gyzr,gzzr, & + avg_detr,densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j,k),tauplus(i,j,k),& + rhoplus(i,j,k),& + velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),& + epsplus(i,j,k),pressplus(i,j,k),Bvecxplus(i,j,k), & + Bvecyplus(i,j,k),Bveczplus(i,j,k),w_lorentzplus(i, j, k)) + + if (densminus(i,j,k) .ne. densminus(i,j,k)) then + !$OMP CRITICAL + call CCTK_WARN(1, "NaN in densminus(i,j,k) (Prim2Con)") + !$OMP END CRITICAL + endif + end do + end do + end do + !$OMP END PARALLEL DO +end subroutine Prim2ConservativePolytypeAM + + /*@@ + @routine prim2conpolytypeAM + @date Aug 31, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke + @desc + Converts from primitive to conservative at a single point + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine prim2conpolytypeAM(handle, gxx, gxy, gxz, gyy, gyz, & + gzz, det, ddens, dsx, dsy, dsz, dtau, & + drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, w) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, & + drho, dvelx, dvely, dvelz, deps, dpress, dBvcx, dBvcy, dBvcz, & + w_tmp, w, vlowx, vlowy, vlowz, sqrtdet + CCTK_INT :: handle + CCTK_REAL :: Bdotv,ab0,b2,blowx,blowy,blowz + character(len=256) NaN_WarnLine + +! begin EOS Omni vars + integer :: n, keytemp, anyerr, keyerr(1) + real*8 :: xpress,xeps,xtemp,xye + n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 + xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 +! end EOS Omni vars + + w_tmp = DOT2(dvelx,dvely,dvelz) + + if (w_tmp .ge. 1.d0) then + ! In theory this should not happen, and even when accepting the fact + ! that numerically it can, one might be tempted to set w to some large + ! value in that case. However, this would lead to completely bogus + ! and hard to trace wrong values below. There is no good value to + ! choose in this case, but something small is probably the best of + ! all bad choices. + !$OMP CRITICAL + write(NaN_WarnLine,'(a80,2g15.6)') 'Infinite Lorentz factor reset. rho, w_tmp: ', drho, w_tmp + call CCTK_WARN(GRHydro_NaN_verbose, NaN_WarnLine) + !$OMP END CRITICAL + w = 1.d-20 + else + w = 1.d0 / sqrt(1.d0 - w_tmp) + endif + + call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,xeps,xtemp,xye,dpress,keyerr,anyerr) + + call EOS_Omni_EpsFromPress(handle,keytemp,GRHydro_eos_rf_prec,n,& + drho,xeps,xtemp,xye,dpress,deps,keyerr,anyerr) + + vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz + vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz + vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz + + Bdotv=DOT(dvelx,dvely,dvelz,dBvcx,dBvcy,dBvcz) + ab0=w*Bdotv + b2 = DOT2(dBvcx,dBvcy,dBvcz)/w**2+Bdotv**2 + blowx = (gxx*dBvcx + gxy*dBvcy + gxz*dBvcz)/w + & + w*Bdotv*vlowx + blowy = (gxy*dBvcx + gyy*dBvcy + gyz*dBvcz)/w + & + w*Bdotv*vlowy + blowz = (gxz*dBvcx + gyz*dBvcy + gzz*dBvcz)/w + & + w*Bdotv*vlowz + + ddens = sqrt(det) * drho * w + dsx = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowx - & + ab0*blowx) + dsy = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowy - & + ab0*blowy) + dsz = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w * vlowz - & + ab0*blowz) + dtau = sqrt(det) * ((drho*(1+deps)+dpress+b2)*w*w - dpress-b2/2.0-ab0**2) - ddens + +end subroutine prim2conpolytypeAM + + + /*@@ + @routine Primitive2ConservativePolyCellsAM + @date Aug 31, 2010 + @author + @desc + Wrapper function that converts primitive to conservative at the + cell centres. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + +subroutine Primitive2ConservativePolyCellsAM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i, j, k + CCTK_REAL :: det + + !$OMP PARALLEL DO PRIVATE(k,j,i,det) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call prim2conpolytypeAM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),& + rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),Bvecx(i,j,k), Bvecy(i,j,k), & + Bvecz(i,j,k), w_lorentz(i,j,k)) + + end do + end do + end do + !$OMP END PARALLEL DO + + if(evolve_Y_e.ne.0) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 + do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 + do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 + Y_e_con(i,j,k) = Y_e(i,j,k) * dens(i,j,k) + enddo + enddo + enddo + !$OMP END PARALLEL DO + endif + + +end subroutine Primitive2ConservativePolyCellsAM + +!!$ +!!$ Prim2Con doesn't change for tracers with the addition of a B-field! +!!$ + diff --git a/src/GRHydro_RegisterVars.cc b/src/GRHydro_RegisterVars.cc index 95a7053..55070ee 100644 --- a/src/GRHydro_RegisterVars.cc +++ b/src/GRHydro_RegisterVars.cc @@ -94,6 +94,12 @@ extern "C" void GRHydro_Register(CCTK_ARGUMENTS) if(clean_divergence) { register_evolved("GRHydro::psidc" , "GRHydro::psidcrhs"); } + } else if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec")) { + register_constrained("HydroBase::Bvec"); + register_evolved("HydroBase::Avec", "GRHydro::Avecrhs"); + if ( CCTK_EQUALS(Avec_gauge, "lorenz")) { + register_evolved("HydroBase::Aphi", "GRHydro::Aphirhs"); + } } // entropycons if(CCTK_EQUALS(entropy_evolution_method,"GRHydro")){ diff --git a/src/GRHydro_RiemannSolveAM.F90 b/src/GRHydro_RiemannSolveAM.F90 new file mode 100644 index 0000000..b51105c --- /dev/null +++ b/src/GRHydro_RiemannSolveAM.F90 @@ -0,0 +1,144 @@ + /*@@ + @file GRHydro_RiemannSolveAM.F90 + @date Sep 1, 2010 + @author + @desc + A wrapper routine to call the correct Riemann solver + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + + /*@@ + @routine RiemannSolveAM + @date Sep 1, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke + @desc + A wrapper routine to switch between the different Riemann solvers. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine RiemannSolveAM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i,j,k + + if (CCTK_EQUALS(riemann_solver,"HLLE").or.CCTK_EQUALS(riemann_solver,"LLF")) then + + call GRHydro_HLLE_AM(CCTK_PASS_FTOF) + + if (evolve_tracer .ne. 0) then + +!!$ There are no special calls for tracers, which care not one whit about B-fields! +!!$ Just call the standard version... + + call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) + + end if + +!!$ else if (CCTK_EQUALS(riemann_solver,"Roe")) then +!!$ +!!$ call GRHydro_RoeSolveAM(CCTK_PASS_FTOF) +!!$ +!!$ if (evolve_tracer .ne. 0) then +!!$ +!!$ call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) +!!$ +!!$ end if +!!$ +!!$ else if (CCTK_EQUALS(riemann_solver,"Marquina")) then +!!$ +!!$ call GRHydro_MarquinaM(CCTK_PASS_FTOF) + +!!$ Tracers are built directly in to the Marquina solver + + else + + call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!") + + end if + +end subroutine RiemannSolveAM + + /*@@ + @routine RiemannSolvePolytypeAM + @date Sep 1, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + The same as above, just specializing to polytropic type EOS. + Currently there is no point to this routine right now. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + + +subroutine RiemannSolvePolytypeAM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i,j,k + + if (CCTK_EQUALS(riemann_solver,"HLLE")) then + + call GRHydro_HLLE_AM(CCTK_PASS_FTOF) + + if (evolve_tracer .ne. 0) then + +!!$ Call the non-MHD version - see above + + call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) + + end if + +!!$ else if (CCTK_EQUALS(riemann_solver,"Roe")) then +!!$ +!!$ call GRHydro_RoeSolve(CCTK_PASS_FTOF) +!!$ +!!$ if (evolve_tracer .ne. 0) then +!!$ +!!$ call GRHydro_HLLE_Tracer(CCTK_PASS_FTOF) +!!$ +!!$ end if +!!$ +!!$ else if (CCTK_EQUALS(riemann_solver,"Marquina")) then +!!$ +!!$ call GRHydro_Marquina(CCTK_PASS_FTOF) + +!!$ Tracers are built directly in to the Marquina solver + + else + + call CCTK_WARN(0, "Roe and Marquina not implemented in MHD yet!!!") + + end if + +end subroutine RiemannSolvePolytypeAM + + + + + diff --git a/src/GRHydro_SourceAM.F90 b/src/GRHydro_SourceAM.F90 new file mode 100644 index 0000000..1e0b190 --- /dev/null +++ b/src/GRHydro_SourceAM.F90 @@ -0,0 +1,679 @@ + /*@@ + @file GRHydro_SourceM.F90 + @date Oct 11, 2010 + @author Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + The geometric source terms for the matter evolution + @enddesc + @@*/ + +! Second order f.d. + +#define DIFF_X_2(q) (0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx) +#define DIFF_Y_2(q) (0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy) +#define DIFF_Z_2(q) (0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz) + +! Fourth order f.d. + +#define DIFF_X_4(q) ((-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \ + q(i-2,j,k)) / 12.d0 * idx) +#define DIFF_Y_4(q) ((-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \ + q(i,j-2,k)) / 12.d0 * idy) +#define DIFF_Z_4(q) ((-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \ + q(i,j,k-2)) / 12.d0 * idz) + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "GRHydro_Macros.h" + +#define velx(i,j,k) vup(i,j,k,1) +#define vely(i,j,k) vup(i,j,k,2) +#define velz(i,j,k) vup(i,j,k,3) +#define Bvecx(i,j,k) Bprim(i,j,k,1) +#define Bvecy(i,j,k) Bprim(i,j,k,2) +#define Bvecz(i,j,k) Bprim(i,j,k,3) +#define Avecx(i,j,k) Avec(i,j,k,1) +#define Avecy(i,j,k) Avec(i,j,k,2) +#define Avecz(i,j,k) Avec(i,j,k,3) +#define Avecrhsx(i,j,k) Avecrhs(i,j,k,1) +#define Avecrhsy(i,j,k) Avecrhs(i,j,k,2) +#define Avecrhsz(i,j,k) Avecrhs(i,j,k,3) + + /*@@ + @routine SourceTermsAM + @date Aug 30, 2010 + @author Tanja Bode, Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + Calculate the geometric source terms and add to the update GFs + @enddesc + @calls + @calledby + @history + Minor alterations of routine from GR3D. + @endhistory + +@@*/ + +subroutine SourceTermsAM(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET kaa, kab, kac, kbb, kbc, kcc + TARGET kxx, kxy, kxz, kyy, kyz, kzz + TARGET betaa, betab, betac + TARGET betax, betay, betaz + TARGET lvel, vel + TARGET lBvec, Bvec + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: 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 :: sqrtdet, det, uxx, uxy, uxz, uyy, uyz, uzz + 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, 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 :: localgxx,localgxy,localgxz,localgyy,localgyz,localgzz + 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, idx, idy, idz + CCTK_REAL :: shiftshiftk, shiftkx, shiftky, shiftkz + CCTK_REAL :: sumTK + CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz + CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz + CCTK_REAL :: invalp, invalp2 + CCTK_REAL :: Avcx_source, Avcy_source, Avcz_source + CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet + CCTK_REAL :: gdg_x, gdg_y, gdg_z !! g^{ik} d_k g_{ij} + + CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum1,dum2,bxlow,bylow,bzlow + CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar + + logical, allocatable, dimension (:,:,:) :: force_spatial_second_order + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: k11, k12, k13, k22, k23, k33 + CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + k11 => kaa + k12 => kab + k13 => kac + k22 => kbb + k23 => kbc + k33 => kcc + beta1 => betaa + beta2 => betab + beta3 => betac + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + k11 => kxx + k12 => kxy + k13 => kxz + k22 => kyy + k23 => kyz + k33 => kzz + beta1 => betax + beta2 => betay + beta3 => betaz + vup => vel + Bprim => Bvec + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define kxx faulty_kxx +#define kxy faulty_kxy +#define kxz faulty_kxz +#define kyy faulty_kyy +#define kyz faulty_kyz +#define kzz faulty_kzz +#define betax faulty_betax +#define betay faulty_betay +#define betaz faulty_betaz +#define gaa faulty_gaa +#define gab faulty_gab +#define gac faulty_gac +#define gbb faulty_gbb +#define gbc faulty_gbc +#define gcc faulty_gcc +#define kaa faulty_kaa +#define kab faulty_kab +#define kac faulty_kac +#define kbb faulty_kbb +#define kbc faulty_kbc +#define kcc faulty_kcc +#define betaa faulty_betaa +#define betab faulty_betab +#define betac faulty_betac +#define vel faulty_vel +#define Bvec faulty_Bvec + + 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) + idx = 1.d0/dx + idy = 1.d0/dy + idz = 1.d0/dz + +!!$ Initialize the update terms to be zero. +!!$ This will guarantee that no garbage in the boundaries is updated. + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + Avecrhs = 0.d0 + + if (evolve_tracer .ne. 0) then + cons_tracerrhs = 0.d0 + end if + + if (evolve_Y_e .ne. 0) then + Y_e_con_rhs = 0.0d0 + endif + + if (clean_divergence .ne. 0) then + psidcrhs=0.d0 + endif + + if (track_divB .ne. 0) then + divB=0.d0 + endif + + if (transport_constraints .ne. 0) then + Evec = 0.d0 + endif + + +!!$ Set up the array for checking the order. We switch to second order +!!$ differencing at boundaries and near excision regions. +!!$ Copied straight from BSSN. + + allocate (force_spatial_second_order(nx,ny,nz)) + force_spatial_second_order = .FALSE. + + if (spatial_order > 2) then + !$OMP PARALLEL DO PRIVATE(i, j, k) + do k = 1 + GRHydro_stencil, nz - GRHydro_stencil + do j = 1 + GRHydro_stencil, ny - GRHydro_stencil + do i = 1 + GRHydro_stencil, nx - GRHydro_stencil + if ((i < 3).or.(i > cctk_lsh(1) - 2).or. & + (j < 3).or.(j > cctk_lsh(2) - 2).or. & + (k < 3).or.(k > cctk_lsh(3) - 2) ) then + force_spatial_second_order(i,j,k) = .TRUE. + else if ( use_mask > 0 ) then + if (minval(emask(i-2:i+2,j-2:j+2,k-2:k+2)) < 0.75d0) then + force_spatial_second_order(i,j,k) = .TRUE. + end if + end if + end do + end do + end do + !$OMP END PARALLEL DO + end if + + !$OMP PARALLEL DO PRIVATE(i, j, k, local_spatial_order,& + !$OMP localgxx,localgxy,localgxz,localgyy,localgyz,localgzz,& + !$OMP det,sqrtdet,shiftx,shifty,shiftz,& + !$OMP dx_betax,dx_betay,dx_betaz,dy_betax,dy_betay,dy_betaz,& + !$OMP dz_betax,dz_betay,dz_betaz,velxshift,velyshift,velzshift,& + !$OMP vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & + !$OMP bdotv,b2,dum1,dum2,bxlow,bylow,bzlow,bt,bx,by,bz,rhohstarW2,pstar,& + !$OMP t00,t0x,t0y,t0z,txx,txy,txz,tyy,tyz,tzz,& + !$OMP dx_alp,dy_alp,dz_alp,& + !$OMP tau_source,sx_source,sy_source,sz_source,& + !$OMP dx_det_bydet,dy_det_bydet,dz_det_bydet,& + !$OMP gdg_x,gdg_y,gdg_z,& + !$OMP Avcx_source,Avcy_source,Avcz_source,evolve_Lorenz_gge,& + !$OMP uxx, uxy, uxz, uyy, uyz, uzz,& + !$OMP dx_gxx, dx_gxy, dx_gxz, dx_gyy, dx_gyz, dx_gzz,& + !$OMP dy_gxx, dy_gxy, dy_gxz, dy_gyy, dy_gyz, dy_gzz,& + !$OMP dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz,& + !$OMP shiftshiftk,shiftkx,shiftky,shiftkz,& + !$OMP sumTK,halfshiftdgx,halfshiftdgy,halfshiftdgz,& + !$OMP halfTdgx,halfTdgy,halfTdgz,invalp,invalp2) + do k=1 + GRHydro_stencil,nz - GRHydro_stencil + do j=1 + GRHydro_stencil,ny - GRHydro_stencil + do i=1 + GRHydro_stencil,nx - GRHydro_stencil + + local_spatial_order = spatial_order + if (force_spatial_second_order(i,j,k)) then + local_spatial_order = 2 + end if + +!!$ Set the metric terms. + + localgxx = g11(i,j,k) + localgxy = g12(i,j,k) + localgxz = g13(i,j,k) + localgyy = g22(i,j,k) + localgyz = g23(i,j,k) + localgzz = g33(i,j,k) + + det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\ + localgyy, localgyz, localgzz) + sqrtdet = sqrt(det) + call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,& + localgxy, localgxz, localgyy, localgyz, localgzz) + + + shiftx = beta1(i,j,k) + shifty = beta2(i,j,k) + shiftz = beta3(i,j,k) + + if (local_spatial_order .eq. 2) then + + dx_betax = DIFF_X_2(beta1) + dx_betay = DIFF_X_2(beta2) + dx_betaz = DIFF_X_2(beta3) + + dy_betax = DIFF_Y_2(beta1) + dy_betay = DIFF_Y_2(beta2) + dy_betaz = DIFF_Y_2(beta3) + + dz_betax = DIFF_Z_2(beta1) + dz_betay = DIFF_Z_2(beta2) + dz_betaz = DIFF_Z_2(beta3) + + else + + dx_betax = DIFF_X_4(beta1) + dx_betay = DIFF_X_4(beta2) + dx_betaz = DIFF_X_4(beta3) + + dy_betax = DIFF_Y_4(beta1) + dy_betay = DIFF_Y_4(beta2) + dy_betaz = DIFF_Y_4(beta3) + + dz_betax = DIFF_Z_4(beta1) + dz_betay = DIFF_Z_4(beta2) + dz_betaz = DIFF_Z_4(beta3) + + end if + + invalp = 1.0d0 / alp(i,j,k) + invalp2 = invalp**2 + velxshift = velx(i,j,k) - shiftx*invalp + velyshift = vely(i,j,k) - shifty*invalp + velzshift = velz(i,j,k) - shiftz*invalp + + call calc_vlow_blow(localgxx,localgxy,localgxz,localgyy,localgyz,localgzz, & + velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k), & + vlowx,vlowy,vlowz,Bvecxlow,Bvecylow,Bveczlow, & + bdotv,b2,dum1,dum2,bxlow,bylow,bzlow) + +!!$ These are the contravariant components + bt = w_lorentz(i,j,k)/alp(i,j,k)*bdotv + bx = Bvecx(i,j,k)/w_lorentz(i,j,k)+w_lorentz(i,j,k)*bdotv*velxshift + by = Bvecy(i,j,k)/w_lorentz(i,j,k)+w_lorentz(i,j,k)*bdotv*velyshift + bz = Bvecz(i,j,k)/w_lorentz(i,j,k)+w_lorentz(i,j,k)*bdotv*velzshift + + rhohstarW2 = (rho(i,j,k)*(one + eps(i,j,k)) + press(i,j,k)+ b2)*& + w_lorentz(i,j,k)**2 + pstar = press(i,j,k)+0.5d0*b2 + +!!$ For a change, these are T^{ij} + + t00 = (rhohstarW2 - pstar)*invalp2-bt**2 + t0x = rhohstarW2*velxshift*invalp +& + pstar*shiftx*invalp2-bt*bx + t0y = rhohstarW2*velyshift*invalp +& + pstar*shifty*invalp2-bt*by + t0z = rhohstarW2*velzshift*invalp +& + pstar*shiftz*invalp2-bt*bz + txx = rhohstarW2*velxshift*velxshift +& + pstar*(uxx - shiftx*shiftx*invalp2)-bx**2 + txy = rhohstarW2*velxshift*velyshift +& + pstar*(uxy - shiftx*shifty*invalp2)-bx*by + txz = rhohstarW2*velxshift*velzshift +& + pstar*(uxz - shiftx*shiftz*invalp2)-bx*bz + tyy = rhohstarW2*velyshift*velyshift +& + pstar*(uyy - shifty*shifty*invalp2)-by**2 + tyz = rhohstarW2*velyshift*velzshift +& + pstar*(uyz - shifty*shiftz*invalp2)-by*bz + tzz = rhohstarW2*velzshift*velzshift +& + pstar*(uzz - shiftz*shiftz*invalp2)-bz**2 + +!!$ Derivatives of the lapse, metric and shift + + if (local_spatial_order .eq. 2) then + + dx_alp = DIFF_X_2(alp) + dy_alp = DIFF_Y_2(alp) + dz_alp = DIFF_Z_2(alp) + + else + + dx_alp = DIFF_X_4(alp) + dy_alp = DIFF_Y_4(alp) + dz_alp = DIFF_Z_4(alp) + + end if + + if (local_spatial_order .eq. 2) then + + dx_gxx = DIFF_X_2(g11) + dx_gxy = DIFF_X_2(g12) + dx_gxz = DIFF_X_2(g13) + dx_gyy = DIFF_X_2(g22) + dx_gyz = DIFF_X_2(g23) + dx_gzz = DIFF_X_2(g33) + dy_gxx = DIFF_Y_2(g11) + dy_gxy = DIFF_Y_2(g12) + dy_gxz = DIFF_Y_2(g13) + dy_gyy = DIFF_Y_2(g22) + dy_gyz = DIFF_Y_2(g23) + dy_gzz = DIFF_Y_2(g33) + dz_gxx = DIFF_Z_2(g11) + dz_gxy = DIFF_Z_2(g12) + dz_gxz = DIFF_Z_2(g13) + dz_gyy = DIFF_Z_2(g22) + dz_gyz = DIFF_Z_2(g23) + dz_gzz = DIFF_Z_2(g33) + + else + + dx_gxx = DIFF_X_4(g11) + dx_gxy = DIFF_X_4(g12) + dx_gxz = DIFF_X_4(g13) + dx_gyy = DIFF_X_4(g22) + dx_gyz = DIFF_X_4(g23) + dx_gzz = DIFF_X_4(g33) + dy_gxx = DIFF_Y_4(g11) + dy_gxy = DIFF_Y_4(g12) + dy_gxz = DIFF_Y_4(g13) + dy_gyy = DIFF_Y_4(g22) + dy_gyz = DIFF_Y_4(g23) + dy_gzz = DIFF_Y_4(g33) + dz_gxx = DIFF_Z_4(g11) + dz_gxy = DIFF_Z_4(g12) + dz_gxz = DIFF_Z_4(g13) + dz_gyy = DIFF_Z_4(g22) + dz_gyz = DIFF_Z_4(g23) + dz_gzz = DIFF_Z_4(g33) + + end if + +!!$ Contract the shift with the extrinsic curvature + + shiftshiftk = shiftx*shiftx*k11(i,j,k) + & + shifty*shifty*k22(i,j,k) + & + shiftz*shiftz*k33(i,j,k) + & + two*(shiftx*shifty*k12(i,j,k) + & + shiftx*shiftz*k13(i,j,k) + & + shifty*shiftz*k23(i,j,k)) + + shiftkx = shiftx*k11(i,j,k) + shifty*k12(i,j,k) + shiftz*k13(i,j,k) + shiftky = shiftx*k12(i,j,k) + shifty*k22(i,j,k) + shiftz*k23(i,j,k) + shiftkz = shiftx*k13(i,j,k) + shifty*k23(i,j,k) + shiftz*k33(i,j,k) + +!!$ Contract the matter terms with the extrinsic curvature + + sumTK = txx*k11(i,j,k) + tyy*k22(i,j,k) + tzz*k33(i,j,k) & + + two*(txy*k12(i,j,k) + txz*k13(i,j,k) + tyz*k23(i,j,k)) + +!!$ Update term for tau + + tau_source = t00* & + (shiftshiftk - (shiftx*dx_alp + shifty*dy_alp + shiftz*dz_alp) )& + + t0x*(-dx_alp + two*shiftkx) & + + t0y*(-dy_alp + two*shiftky) & + + t0z*(-dz_alp + two*shiftkz) & + + sumTK + +!!$ The following looks very little like the terms in the +!!$ standard papers. Take a look in the ThornGuide to see why +!!$ it is really the same thing. + +!!$ Contract the shift with derivatives of the metric + + halfshiftdgx = half*(shiftx*shiftx*dx_gxx + & + shifty*shifty*dx_gyy + shiftz*shiftz*dx_gzz) + & + shiftx*shifty*dx_gxy + shiftx*shiftz*dx_gxz + & + shifty*shiftz*dx_gyz + halfshiftdgy = half*(shiftx*shiftx*dy_gxx + & + shifty*shifty*dy_gyy + shiftz*shiftz*dy_gzz) + & + shiftx*shifty*dy_gxy + shiftx*shiftz*dy_gxz + & + shifty*shiftz*dy_gyz + halfshiftdgz = half*(shiftx*shiftx*dz_gxx + & + shifty*shifty*dz_gyy + shiftz*shiftz*dz_gzz) + & + shiftx*shifty*dz_gxy + shiftx*shiftz*dz_gxz + & + shifty*shiftz*dz_gyz + +!!$ Contract the matter with derivatives of the metric + + halfTdgx = half*(txx*dx_gxx + tyy*dx_gyy + tzz*dx_gzz) +& + txy*dx_gxy + txz*dx_gxz + tyz*dx_gyz + halfTdgy = half*(txx*dy_gxx + tyy*dy_gyy + tzz*dy_gzz) +& + txy*dy_gxy + txz*dy_gxz + tyz*dy_gyz + halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) +& + txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz + + + + + sx_source = t00*& + (halfshiftdgx - alp(i,j,k)*dx_alp) + halfTdgx + & + 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) +& + rhohstarW2*invalp*(vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz) -& + bt*(bxlow*dx_betax + bylow*dx_betay + bzlow*dx_betaz) + + sy_source = t00*& + (halfshiftdgy - alp(i,j,k)*dy_alp) + halfTdgy + & + 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) +& + rhohstarW2*invalp*(vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz) -& + bt*(bxlow*dy_betax + bylow*dy_betay + bzlow*dy_betaz) + + sz_source = t00*& + (halfshiftdgz - alp(i,j,k)*dz_alp) + halfTdgz + & + 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) +& + rhohstarW2*invalp*(vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz) -& + bt*(bxlow*dz_betax + bylow*dz_betay + bzlow*dz_betaz) + + !! B^i and A^i both live in cell centers currently + Avcx_source = sqrtdet*(vely(i,j,k)*Bvecz(i,j,k) - velz(i,j,k)*Bvecy(i,j,k)) + Avcy_source = sqrtdet*(velz(i,j,k)*Bvecx(i,j,k) - velx(i,j,k)*Bvecz(i,j,k)) + Avcz_source = sqrtdet*(velx(i,j,k)*Bvecy(i,j,k) - vely(i,j,k)*Bvecx(i,j,k)) + + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = 0.d0 + end if + + densrhs(i,j,k) = 0.d0 + srhs(i,j,k,1) = alp(i,j,k)*sqrtdet*sx_source + srhs(i,j,k,2) = alp(i,j,k)*sqrtdet*sy_source + srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source + taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source + Avecrhsx(i,j,k) = Avcx_source + Avecrhsy(i,j,k) = Avcy_source + Avecrhsz(i,j,k) = Avcz_source + + enddo + enddo + enddo + !$OMP END PARALLEL DO + + deallocate(force_spatial_second_order) + +#if(0) /* poison edges of domain */ + if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then + last_iteration_seen = cctk_iteration + reflevel = grhydro_reflevel + mol_substep = 0 + else + mol_substep = mol_substep + 1 + end if + do k = 1, GRHydro_stencil*mol_substep + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Avecx(i,j,k) = -1d100 + Avecy(i,j,k) = -1d100 + Avecz(i,j,k) = -1d100 + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = -1d100 + end if + end do + end do + end do + do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Avecx(i,j,k) = -1d100 + Avecy(i,j,k) = -1d100 + Avecz(i,j,k) = -1d100 + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = -1d100 + end if + end do + end do + end do + do i = 1, GRHydro_stencil*mol_substep + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Avecx(i,j,k) = -1d100 + Avecy(i,j,k) = -1d100 + Avecz(i,j,k) = -1d100 + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = -1d100 + end if + end do + end do + end do + do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Avecx(i,j,k) = -1d100 + Avecy(i,j,k) = -1d100 + Avecz(i,j,k) = -1d100 + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = -1d100 + end if + end do + end do + end do + do j = 1, GRHydro_stencil*mol_substep + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Avecx(i,j,k) = -1d100 + Avecy(i,j,k) = -1d100 + Avecz(i,j,k) = -1d100 + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = -1d100 + end if + end do + end do + end do + do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + do k = 1, cctk_lsh(3) + dens(i,j,k) = -1d100 + Scon(i,j,k,1) = -1d100 + Scon(i,j,k,2) = -1d100 + Scon(i,j,k,3) = -1d100 + tau(i,j,k) = -1d100 + Avecx(i,j,k) = -1d100 + Avecy(i,j,k) = -1d100 + Avecz(i,j,k) = -1d100 + if ( evolve_Lorenz_gge.gt.0 ) then + Aphi(i,j,k) = -1d100 + end if + end do + end do + end do +#endif + +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_vel +#undef faulty_Bvec +#undef faulty_gxx_p +#undef faulty_gxy_p +#undef faulty_gxz_p +#undef faulty_gyy_p +#undef faulty_gyz_p +#undef faulty_gzz_p +#undef faulty_vel_p +#undef faulty_Bvec_p +#undef faulty_gxx_p_p +#undef faulty_gxy_p_p +#undef faulty_gxz_p_p +#undef faulty_gyy_p_p +#undef faulty_gyz_p_p +#undef faulty_gzz_p_p +#undef faulty_vel_p_p +#undef faulty_Bvec_p_p +end subroutine SourceTermsAM + + + diff --git a/src/GRHydro_UpdateMaskM.F90 b/src/GRHydro_UpdateMaskM.F90 index a965f53..84fdec5 100644 --- a/src/GRHydro_UpdateMaskM.F90 +++ b/src/GRHydro_UpdateMaskM.F90 @@ -503,3 +503,478 @@ subroutine GRHydro_InitialAtmosphereResetM(CCTK_ARGUMENTS) end subroutine GRHydro_InitialAtmosphereResetM + + /*@@ + @routine GRHydro_AtmosphereResetAM + @date Sep 2, 2010 + @author Tanja Bode, Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + After MoL has evolved, if a point is supposed to be reset then do so. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +subroutine GRHydro_AtmosphereResetAM(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i, j, k + CCTK_REAL :: det + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,xeps,xtemp,xye +! end EOS Omni vars + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + +! begin EOS Omni vars + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0 +! end EOS Omni vars + + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define vel faulty_vel +#define Bvec faulty_Bvec + + if (verbose.eq.1) call CCTK_INFO("Entering AtmosphereReset.") + +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + + if (atmosphere_mask(i, j, k) .ne. 0) then + + rho(i,j,k) = GRHydro_rho_min + vup(i,j,k,1) = 0.0d0 + vup(i,j,k,2) = 0.0d0 + vup(i,j,k,3) = 0.0d0 + det = SPATIAL_DETERMINANT(g11(i,j,k), g12(i,j,k), g13(i,j,k), \ + g22(i,j,k), g23(i,j,k), g33(i,j,k)) + + if(evolve_temper.ne.0) then + + ! set the temperature to be relatively low + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),& + g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),& + det, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3),& + tau(i,j,k),Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),& + eps(i,j,k),press(i,j,k),Bprim(i,j,k,1), & + Bprim(i,j,k,2), Bprim(i,j,k,3), w_lorentz(i,j,k),& + temperature(i,j,k),y_e(i,j,k)) + y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k) + + else + + call prim2conpolytypeM(GRHydro_polytrope_handle, & + g11(i,j,k), g12(i,j,k), g13(i,j,k), & + g22(i,j,k), g23(i,j,k), g33(i,j,k), det, & + dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + tau(i,j,k), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k), vup(i,j,k,1), vup(i,j,k,2), & + vup(i,j,k,3), eps(i,j,k), press(i,j,k), & + Bprim(i,j,k,1),Bprim(i,j,k,2),Bprim(i,j,k,3),w_lorentz(i,j,k)) + if (wk_atmosphere .eq. 0) then + atmosphere_mask(i, j, k) = 0 + atmosphere_mask_real(i, j, k) = 0 + end if + + end if + endif + + end do + end do + end do +!$OMP END PARALLEL DO + +!!$ call GRHydro_BoundariesM(CCTK_PASS_FTOF) +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_vel +#undef faulty_Bvec + +end subroutine GRHydro_AtmosphereResetAM + +subroutine GRHydro_InitialAtmosphereResetAM(CCTK_ARGUMENTS) + + implicit none + + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET gaa, gab, gac, gbb, gbc, gcc + TARGET gxx, gxy, gxz, gyy, gyz, gzz + TARGET lvel, vel + TARGET lBvec, Bvec + TARGET gaa_p, gab_p, gac_p, gbb_p, gbc_p, gcc_p + TARGET gxx_p, gxy_p, gxz_p, gyy_p, gyz_p, gzz_p + TARGET lvel_p, vel_p + TARGET lBvec_p, Bvec_p + TARGET gaa_p_p, gab_p_p, gac_p_p, gbb_p_p, gbc_p_p, gcc_p_p + TARGET gxx_p_p, gxy_p_p, gxz_p_p, gyy_p_p, gyz_p_p, gzz_p_p + TARGET lvel_p_p, vel_p_p + TARGET lBvec_p_p, Bvec_p_p + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + CCTK_INT :: i, j, k + CCTK_REAL :: det + CCTK_REAL :: rho_min + + CCTK_INT :: eos_handle + + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup, Bprim + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11_p, g12_p, g13_p, g22_p, & + g23_p, g33_p + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup_p, Bprim_p + CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11_p_p, g12_p_p, g13_p_p, g22_p_p, & + g23_p_p, g33_p_p + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup_p_p, Bprim_p_p + +! begin EOS Omni vars + integer :: n,keytemp,anyerr,keyerr(1) + real*8 :: xpress,xeps,xtemp,xye + n=1;keytemp=0;anyerr=0;keyerr(1)=0 + xpress=0.0d0;xye=0.0d0;xeps=0.0d0;xtemp=0.0d0 +! end EOS Omni vars + + ! save memory when MP is not used + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11 => gaa + g12 => gab + g13 => gac + g22 => gbb + g23 => gbc + g33 => gcc + vup => lvel + Bprim => lBvec + else + g11 => gxx + g12 => gxy + g13 => gxz + g22 => gyy + g23 => gyz + g33 => gzz + vup => vel + Bprim => Bvec + end if + if (timelevels .gt. 1) then + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11_p => gaa_p + g12_p => gab_p + g13_p => gac_p + g22_p => gbb_p + g23_p => gbc_p + g33_p => gcc_p + vup_p => lvel_p + Bprim_p => Bvec_p + else + g11_p => gxx_p + g12_p => gxy_p + g13_p => gxz_p + g22_p => gyy_p + g23_p => gyz_p + g33_p => gzz_p + vup_p => vel_p + Bprim_p => Bvec_p + end if + end if + if (timelevels .gt. 2) then + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + g11_p_p => gaa_p_p + g12_p_p => gab_p_p + g13_p_p => gac_p_p + g22_p_p => gbb_p_p + g23_p_p => gbc_p_p + g33_p_p => gcc_p_p + vup_p_p => lvel_p_p + Bprim_p_p => lBvec_p_p + else + g11_p_p => gxx_p_p + g12_p_p => gxy_p_p + g13_p_p => gxz_p_p + g22_p_p => gyy_p_p + g23_p_p => gyz_p_p + g33_p_p => gzz_p_p + vup_p_p => vel_p_p + Bprim_p_p => Bvec_p_p + end if + end if +#define gxx faulty_gxx +#define gxy faulty_gxy +#define gxz faulty_gxz +#define gyy faulty_gyy +#define gyz faulty_gyz +#define gzz faulty_gzz +#define vel faulty_vel +#define Bvec faulty_Bvec +#define gxx_p faulty_gxx_p +#define gxy_p faulty_gxy_p +#define gxz_p faulty_gxz_p +#define gyy_p faulty_gyy_p +#define gyz_p faulty_gyz_p +#define gzz_p faulty_gzz_p +#define vel_p faulty_vel_p +#define Bvec_p faulty_Bvec_p +#define gxx_p_p faulty_gxx_p_p +#define gxy_p_p faulty_gxy_p_p +#define gxz_p_p faulty_gxz_p_p +#define gyy_p_p faulty_gyy_p_p +#define gyz_p_p faulty_gyz_p_p +#define gzz_p_p faulty_gzz_p_p +#define vel_p_p faulty_vel_p_p +#define Bvec_p_p faulty_Bvec_p_p + + + eos_handle = GRHydro_polytrope_handle + + rho_min = GRHydro_rho_min + if (initial_atmosphere_factor .gt. 0) then + rho_min = rho_min * initial_atmosphere_factor + endif + +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) + do k = 1, cctk_lsh(3) + do j = 1, cctk_lsh(2) + do i = 1, cctk_lsh(1) + + if (rho(i,j,k) .le. rho_min .or. & + (GRHydro_enable_internal_excision .ne. 0 .and. & + hydro_excision_mask(i,j,k) .ne. 0) ) then + rho(i,j,k) = rho_min + vup(i,j,k,1) = 0.0d0 + vup(i,j,k,2) = 0.0d0 + vup(i,j,k,3) = 0.0d0 + + det = SPATIAL_DETERMINANT(g11(i,j,k), g12(i,j,k), g13(i,j,k), \ + g22(i,j,k), g23(i,j,k), g33(i,j,k)) + + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + temperature(i,j,k) = grhydro_hot_atmo_temp + y_e(i,j,k) = grhydro_hot_atmo_Y_e + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),& + press(i,j,k),keyerr,anyerr) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11(i,j,k),& + g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),& + det, dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3),& + tau(i,j,k),Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),& + eps(i,j,k),press(i,j,k),Bprim(i,j,k,1), & + Bprim(i,j,k,2), Bprim(i,j,k,3), w_lorentz(i,j,k),& + temperature(i,j,k),y_e(i,j,k)) + y_e_con(i,j,k) = dens(i,j,k) * y_e(i,j,k) + + else + + keytemp = 0 + call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),xeps,xtemp,xye,press(i,j,k),eps(i,j,k),keyerr,anyerr) + + call prim2conpolytypeM(eos_handle, & + g11(i,j,k), g12(i,j,k), g13(i,j,k), & + g22(i,j,k), g23(i,j,k), g33(i,j,k), det, & + dens(i,j,k), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + tau(i,j,k), Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),& + rho(i,j,k), vup(i,j,k,1), vup(i,j,k,2), & + vup(i,j,k,3), eps(i,j,k), press(i,j,k), & + Bprim(i,j,k,1),Bprim(i,j,k,2),Bprim(i,j,k,3),w_lorentz(i,j,k)) + end if + end if + if (timelevels .gt. 1) then + if (rho_p(i,j,k) .le. rho_min) then + rho_p(i,j,k) = rho_min + vup_p(i,j,k,1) = 0.0d0 + vup_p(i,j,k,2) = 0.0d0 + vup_p(i,j,k,3) = 0.0d0 + + det = SPATIAL_DETERMINANT(g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), \ + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k)) + + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + temperature_p(i,j,k) = grhydro_hot_atmo_temp + y_e_p(i,j,k) = grhydro_hot_atmo_Y_e + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),eps_p(i,j,k),temperature_p(i,j,k),y_e_p(i,j,k),& + press_p(i,j,k),keyerr,anyerr) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p(i,j,k),& + g12_p(i,j,k),g13_p(i,j,k),g22_p(i,j,k),g23_p(i,j,k),g33_p(i,j,k),& + det, dens_p(i,j,k),scon_p(i,j,k,1),scon_p(i,j,k,2),scon_p(i,j,k,3),& + tau_p(i,j,k),Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& + rho_p(i,j,k),vup_p(i,j,k,1),vup_p(i,j,k,2),vup_p(i,j,k,3),& + eps_p(i,j,k),press_p(i,j,k),Bprim_p(i,j,k,1), & + Bprim_p(i,j,k,2), Bprim_p(i,j,k,3), w_lorentz_p(i,j,k),& + temperature_p(i,j,k),y_e_p(i,j,k)) + y_e_con_p(i,j,k) = dens_p(i,j,k) * y_e_p(i,j,k) + + else + + keytemp = 0 + call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),eps_p(i,j,k),xtemp,xye,press_p(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p(i,j,k),xeps,xtemp,xye,press_p(i,j,k),eps_p(i,j,k),keyerr,anyerr) + + call prim2conpolytypeM(eos_handle, & + g11_p(i,j,k), g12_p(i,j,k), g13_p(i,j,k), & + g22_p(i,j,k), g23_p(i,j,k), g33_p(i,j,k), det, & + dens_p(i,j,k), scon_p(i,j,k,1), scon_p(i,j,k,2), scon_p(i,j,k,3), & + tau_p(i,j,k), Bcons_p(i,j,k,1),Bcons_p(i,j,k,2),Bcons_p(i,j,k,3),& + rho_p(i,j,k), vup_p(i,j,k,1), vup_p(i,j,k,2), & + vup_p(i,j,k,3), eps_p(i,j,k), press_p(i,j,k), & + Bprim_p(i,j,k,1),Bprim_p(i,j,k,2),Bprim_p(i,j,k,3),w_lorentz_p(i,j,k)) + end if + end if + end if + + if (timelevels .gt. 2) then + if (rho_p_p(i,j,k) .le. rho_min) then + rho_p_p(i,j,k) = rho_min + vup_p_p(i,j,k,1) = 0.0d0 + vup_p_p(i,j,k,2) = 0.0d0 + vup_p_p(i,j,k,3) = 0.0d0 + + det = SPATIAL_DETERMINANT(g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), \ + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k)) + + if(evolve_temper.ne.0) then +! ! set the temperature to be relatively low + temperature_p_p(i,j,k) = grhydro_hot_atmo_temp + y_e_p_p(i,j,k) = grhydro_hot_atmo_Y_e + keytemp = 1 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),eps_p_p(i,j,k),temperature_p_p(i,j,k),y_e_p_p(i,j,k),& + press_p_p(i,j,k),keyerr,anyerr) + call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),g11_p_p(i,j,k),& + g12_p_p(i,j,k),g13_p_p(i,j,k),g22_p_p(i,j,k),g23_p_p(i,j,k),g33_p_p(i,j,k),& + det, dens_p_p(i,j,k),scon_p_p(i,j,k,1),scon_p_p(i,j,k,2),scon_p_p(i,j,k,3),& + tau_p_p(i,j,k),Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& + rho_p_p(i,j,k),vup_p_p(i,j,k,1),vup_p_p(i,j,k,2),vup_p_p(i,j,k,3),& + eps_p_p(i,j,k),press_p_p(i,j,k),Bprim_p_p(i,j,k,1), & + Bprim_p_p(i,j,k,2), Bprim_p_p(i,j,k,3), w_lorentz_p_p(i,j,k),& + temperature_p_p(i,j,k),y_e_p_p(i,j,k)) + y_e_con_p_p(i,j,k) = dens_p_p(i,j,k) * y_e_p_p(i,j,k) + + else + + keytemp = 0 + call EOS_Omni_press(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),eps_p_p(i,j,k),xtemp,xye,press_p_p(i,j,k),keyerr,anyerr) + call EOS_Omni_EpsFromPress(eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho_p_p(i,j,k),xeps,xtemp,xye,press_p_p(i,j,k),eps_p_p(i,j,k),keyerr,anyerr) + + call prim2conpolytypeM(eos_handle, & + g11_p_p(i,j,k), g12_p_p(i,j,k), g13_p_p(i,j,k), & + g22_p_p(i,j,k), g23_p_p(i,j,k), g33_p_p(i,j,k), det, & + dens_p_p(i,j,k), scon_p_p(i,j,k,1), scon_p_p(i,j,k,2), scon_p_p(i,j,k,3), & + tau_p_p(i,j,k), Bcons_p_p(i,j,k,1),Bcons_p_p(i,j,k,2),Bcons_p_p(i,j,k,3),& + rho_p_p(i,j,k), vup_p_p(i,j,k,1), vup_p_p(i,j,k,2), & + vup_p_p(i,j,k,3), eps_p_p(i,j,k), press_p_p(i,j,k), & + Bprim_p_p(i,j,k,1),Bprim_p_p(i,j,k,2),Bprim_p_p(i,j,k,3),w_lorentz_p_p(i,j,k)) + end if + end if + end if + + end do + end do + end do +!$OMP END PARALLEL DO + + write(*,*) " GRHydro_InitialAtmosphereReset" +!!$ call GRHydro_BoundariesM(CCTK_PASS_FTOF) +#undef faulty_gxx +#undef faulty_gxy +#undef faulty_gxz +#undef faulty_gyy +#undef faulty_gyz +#undef faulty_gzz +#undef faulty_vel +#undef faulty_Bvec +#undef faulty_gxx_p +#undef faulty_gxy_p +#undef faulty_gxz_p +#undef faulty_gyy_p +#undef faulty_gyz_p +#undef faulty_gzz_p +#undef faulty_vel_p +#undef faulty_Bvec_p +#undef faulty_gxx_p_p +#undef faulty_gxy_p_p +#undef faulty_gxz_p_p +#undef faulty_gyy_p_p +#undef faulty_gyz_p_p +#undef faulty_gzz_p_p +#undef faulty_vel_p_p +#undef faulty_Bvec_p_p + +end subroutine GRHydro_InitialAtmosphereResetAM + diff --git a/src/make.code.defn b/src/make.code.defn index 02f56dd..5a1b96e 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -57,6 +57,12 @@ SRCS = Utils.F90 \ GRHydro_TmunuM.F90 \ GRHydro_UpdateMaskM.F90 \ GRHydro_UtilsM.F90 \ + GRHydro_BvecfromAvec.F90 \ + GRHydro_RiemannSolveAM.F90 \ + GRHydro_HLLE_AM.F90 \ + GRHydro_Con2PrimAM.F90 \ + GRHydro_Prim2ConAM.F90 \ + GRHydro_SourceAM.F90 \ GRHydro_TransformTensorBasis.c \ GRHydro_Jacobian_state.c \ GRHydro_PPMMReconstruct_drv.F90\ |