[Commits] [svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 593)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Mon Dec 2 15:50:25 CST 2013
User: rhaas
Date: 2013/12/02 03:50 PM
Modified:
/trunk/src/
GRHydro_Con2PrimM.F90, GRHydro_UpdateMask.F90, GRHydro_UpdateMaskM.F90
Log:
GRHydro: Better handling of excision for MHD.
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Con2PrimM.F90
Delta lines: +8 -43
===================================================================
--- trunk/src/GRHydro_Con2PrimM.F90 2013-11-25 15:39:05 UTC (rev 592)
+++ trunk/src/GRHydro_Con2PrimM.F90 2013-12-02 21:50:25 UTC (rev 593)
@@ -164,10 +164,6 @@
!do not compute if in atmosphere region!
if (atmosphere_mask(i,j,k) .gt. 0) cycle
- !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
sdet = sdetg(i,j,k)
@@ -192,43 +188,6 @@
endif
- ! In excised region, set to atmosphere!
- if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .gt. 0)) then
- SET_ATMO_MIN(dens(i,j,k), sdet*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance)
- SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(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
-
- !!$ 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
-
- 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)
@@ -239,8 +198,8 @@
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
+
+ if ( (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .gt. 0)) .or. (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)
@@ -256,6 +215,10 @@
vup(i,j,k,:) = 0.d0
w_lorentz(i,j,k) = 1.d0
+ Bprim(i,j,k,:) = 0.0d0
+ Bcons(i,j,k,:) = 0.0d0
+ b2 = 0
+
if(evolve_temper.ne.0) then
! set hot atmosphere values
temperature(i,j,k) = grhydro_hot_atmo_temp
@@ -294,6 +257,8 @@
end if
+
+
if(evolve_temper.eq.0) then
File [modified]: GRHydro_UpdateMask.F90
Delta lines: +25 -0
===================================================================
--- trunk/src/GRHydro_UpdateMask.F90 2013-11-25 15:39:05 UTC (rev 592)
+++ trunk/src/GRHydro_UpdateMask.F90 2013-12-02 21:50:25 UTC (rev 593)
@@ -10,6 +10,7 @@
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
#include "GRHydro_Macros.h"
#include "SpaceMask.h"
@@ -30,12 +31,24 @@
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i,j,k
CCTK_REAL :: frac
+ logical :: evolve_Bvec
+ logical :: evolve_Avec
frac = CCTK_DELTA_TIME
+ evolve_Bvec = .false.
+ evolve_Avec = .false.
+
+ if (CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then
+ evolve_Bvec = .true.
+ elseif (CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") .or. CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec_centered")) then
+ evolve_Avec = .true.
+ endif
+
if(evolve_temper.ne.1.and.evolve_Y_e.ne.1) then
!$OMP PARALLEL DO PRIVATE(k,j,i)
do k = 1+GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil
@@ -48,6 +61,12 @@
densrhs(i,j,k) = 0.0d0
srhs(i,j,k,:) = 0.0d0
taurhs(i,j,k) = 0.0d0
+
+ if (evolve_Bvec) then
+ Bconsrhs(i,j,k,:) = 0.0d0
+ endif
+
+ ! TODO: Need to set Avecrhs and Aphirhs to zero for centered / or staggered case if vector potential is evolved
! Set real-valued mask! This will be sync'ed and right after syncing translated to
! our standard integer based mask (so that atmosphere_mask is still valid!).
@@ -70,7 +89,13 @@
densrhs(i,j,k) = 0.0d0
srhs(i,j,k,:) = 0.0d0
taurhs(i,j,k) = 0.0d0
+
+ if (evolve_Bvec) then
+ Bconsrhs(i,j,k,:) = 0.0d0
+ endif
+ ! TODO: Need to set Avecrhs and Aphirhs to zero for centered / or staggered case if vector potential is evolved
+
! Set real-valued mask! This will be sync'ed and right after syncing translated to
! our standard integer based mask (so that atmosphere_mask is still valid!).
atmosphere_mask_real(i,j,k) = 1
File [modified]: GRHydro_UpdateMaskM.F90
Delta lines: +2 -0
===================================================================
--- trunk/src/GRHydro_UpdateMaskM.F90 2013-11-25 15:39:05 UTC (rev 592)
+++ trunk/src/GRHydro_UpdateMaskM.F90 2013-12-02 21:50:25 UTC (rev 593)
@@ -115,6 +115,8 @@
vup(i,j,k,3) = 0.0d0
sdet = sdetg(i,j,k)
+ Bprim(i,j,k,:) = 0.0d0
+
if(evolve_temper.ne.0) then
! set the temperature to be relatively low
More information about the Commits
mailing list