User: rhaas Date: 2013/01/14 08:23 AM

Modified: /trunk/ param.ccl /trunk/src/ GRHydro_Con2Prim.F90, GRHydro_Macros.h, GRHydro_Minima.F90, GRHydro_Reconstruct.F90, GRHydro_UpdateMask.F90

Log: GRHydro: Introduce variable atmopshere level / tolerance as function of radius. Off by default.

From: Christian Reisswig <reisswig@scriwalker.(none)>

File Changes:

Directory: /trunk/src/ ======================

File [modified]: GRHydro_Con2Prim.F90 Delta lines: +28 -20 =================================================================== --- trunk/src/GRHydro_Con2Prim.F90 2013-01-14 14:23:24 UTC (rev 452) +++ trunk/src/GRHydro_Con2Prim.F90 2013-01-14 14:23:27 UTC (rev 453) @@ -56,7 +56,7 @@

integer :: CCTK_MyProc integer :: i, j, k, itracer, nx, ny, nz - CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin + CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin, dummy1, dummy2 logical :: epsnegative character*256 :: warnline

@@ -112,7 +112,7 @@

!$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,& - !$OMP warnline) + !$OMP warnline, dummy1, dummy2) do k = 1, nz do j = 1, ny do i = 1, nx @@ -168,10 +168,11 @@ GRHydro_Y_e_min) endif

- if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + !if ( dens(i,j,k) .le. sqrt(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + IF_BELOW_ATMO(dens(i,j,k), sqrt(det)*GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) then

- dens(i,j,k) = sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) - rho(i,j,k) = GRHydro_rho_min + SET_ATMO_MIN(dens(i,j,k), sqrt(det)*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 @@ -272,8 +273,8 @@ !$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 + dens(i,j,k) = ATMOMIN(sqrt(det)*GRHydro_rho_min, r(i,j,k)) !sqrt(det)*GRHydro_rho_min !/(1.d0+GRHydro_atmo_tolerance) + rho(i,j,k) = ATMOMIN(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 @@ -413,7 +414,7 @@

CCTK_REAL s2, c0, c1, c2, c3, c4, f, df, ftol, v2, w, vlowx, vlowy, vlowz CCTK_INT count, i, handle, GRHydro_reflevel - CCTK_REAL GRHydro_C2P_failed + CCTK_REAL GRHydro_C2P_failed, dummy1, dummy2 CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin character(len=200) warnline @@ -511,7 +512,7 @@

! for safety, let's set the point to atmosphere - rho = GRHydro_rho_min + SET_ATMO_MIN(rho, GRHydro_rho_min, r) udens = rho dens = sqrt(det) * rho pnew = pmin @@ -682,8 +683,9 @@

!!$ Calculate primitive variables from root

- if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then - rho = GRHydro_rho_min + !if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + IF_BELOW_ATMO(rho, GRHydro_rho_min, GRHydro_atmo_tolerance, r) then + SET_ATMO_MIN(rho, GRHydro_rho_min, r) !GRHydro_rho_min udens = rho dens = sqrt(det) * rho ! epsilon = (sqrt( (utau + pnew + udens)**2) - pnew - udens)/udens @@ -777,7 +779,7 @@ CCTK_REAL GRHydro_C2P_failed CCTK_REAL udens, usx, usy, usz, utau, pold, pnew, epsold, epsnew, w2, & w2mhalf, temp1, drhobydpress, depsbydpress, dpressbydeps, dpressbydrho, pmin, epsmin - CCTK_REAL epsminl,pminl,plow,tmp + CCTK_REAL epsminl,pminl,plow,tmp, dummy1, dummy2 CCTK_REAL local_perc_ptol

character(len=256) warnline @@ -1218,8 +1220,9 @@

!!$ Calculate primitive variables from root

- if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then - rho = GRHydro_rho_min + !if (rho .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + IF_BELOW_ATMO(rho, GRHydro_rho_min, GRHydro_atmo_tolerance, r) then + SET_ATMO_MIN(rho, GRHydro_rho_min, r) !GRHydro_rho_min udens = rho dens = sqrt(det) * rho temp = GRHydro_hot_atmo_temp @@ -1660,7 +1663,7 @@ CCTK_REAL s2, f, df, vlowx, vlowy, vlowz CCTK_INT count, handle, GRHydro_reflevel CCTK_REAL udens, usx, usy, usz, rhoold, rhonew, epsold, epsnew, & - enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac, GRHydro_C2P_failed + enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac, GRHydro_C2P_failed, dummy1, dummy2 character(len=200) warnline

! begin EOS Omni vars @@ -1727,7 +1730,7 @@ !$OMP END CRITICAL

! for safety, let's set the point to atmosphere - rhonew = GRHydro_rho_min + SET_ATMO_MIN(rhonew, GRHydro_rho_min, r) udens = rhonew dens = sqrtdet * rhonew

@@ -1792,8 +1795,10 @@

!!$ Calculate primitive variables from root

- if (rhonew .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then - rhonew = GRHydro_rho_min + + !if (rhonew .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then + IF_BELOW_ATMO(rhonew, GRHydro_rho_min, GRHydro_atmo_tolerance, r) then + SET_ATMO_MIN(rhonew, GRHydro_rho_min, r) !GRHydro_rho_min udens = rhonew ! before 2009/10/07 dens was not reset and was used with its negative ! value below; this has since been changed, but the change produces @@ -2232,6 +2237,8 @@

integer :: i, j, k, nx, ny, nz character(len=300) warnline + CCTK_REAL :: dummy1, dummy2 +

! call CCTK_INFO("Checking the C2P failure mask.")

@@ -2242,7 +2249,7 @@ ! do not check if told not to if(GRHydro_reflevel .lt. GRHydro_c2p_warn_from_reflevel) return

- !$OMP PARALLEL DO PRIVATE(i,j,k) + !$OMP PARALLEL DO PRIVATE(i,j,k, dummy1, dummy2) do k = 1, nz do j = 1, ny do i = 1, nx @@ -2256,7 +2263,8 @@ cycle endif

- 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) ) then + IF_BELOW_ATMO(rho(i,j,k), GRHydro_rho_min, GRHydro_atmo_tolerance, r(i,j,k)) then cycle end if

File [modified]: GRHydro_Macros.h Delta lines: +28 -0 =================================================================== --- trunk/src/GRHydro_Macros.h 2013-01-14 14:23:24 UTC (rev 452) +++ trunk/src/GRHydro_Macros.h 2013-01-14 14:23:27 UTC (rev 453) @@ -10,3 +10,31 @@ #define DOTP2(gxx_,gxy_,gxz_,gyy_,gyz_,gzz_,x_,y_,z_) \ ( (gxx_)*(x_)**2+(gyy_)*(y_)**2+(gzz_)*(z_)**2+ \ 2.0*( (gxy_)*(x_)*(y_)+(gxz_)*(x_)*(z_)+(gyz_)*(y_)*(z_) ) ) + + +#define IF_BELOW_ATMO(rho, rho_min, rho_tol, r) \ + dummy1 = atmo_tolerance_radius &&\ + dummy2 = atmo_falloff_radius &&\ + if (r .gt. atmo_tolerance_radius) then &&\ + dummy1 = r &&\ + endif &&\ + if (r .gt. atmo_falloff_radius) then &&\ + dummy2 = r &&\ + endif &&\ + if (rho .le. rho_min*(1.0d0 + rho_tol * (dummy1/atmo_tolerance_radius)**atmo_tolerance_power) * (atmo_falloff_radius/dummy2)**atmo_falloff_power) + + + +!#define ATMOCHECK(rho, rho_min, rho_tol, r) \ +! (atmo_type .eq. 0 .and. rho .le. rho_min*(1.0d0 + rho_tol)) .or. \ +! (atmo_type .eq. 1 .and. rho .le. rho_min*(1.0d0 + rho_tol * (atmo_tolerance_radius/r)**atmo_tolerance_power)) .or. \ +! (atmo_type .eq. 2 .and. rho .le. rho_min*(1.d00 + rho_tol) * (atmo_falloff_radius/r)**atmo_falloff_power) .or. \ +! (atmo_type .eq. 3 .and. rho .le. rho_min*(1.0d0 + rho_tol * (r/atmo_tolerance_radius)**atmo_tolerance_power) * (r/atmo_falloff_radius)**atmo_falloff_power) + + +#define SET_ATMO_MIN(rho_, rho_min, r) &&\ + dummy1 = atmo_falloff_radius &&\ + if (r .gt. atmo_falloff_radius) then &&\ + dummy1 = r &&\ + endif &&\ + rho_ = rho_min * (atmo_falloff_radius/dummy1)**atmo_falloff_power

File [modified]: GRHydro_Minima.F90 Delta lines: +6 -2 =================================================================== --- trunk/src/GRHydro_Minima.F90 2013-01-14 14:23:24 UTC (rev 452) +++ trunk/src/GRHydro_Minima.F90 2013-01-14 14:23:27 UTC (rev 453) @@ -11,6 +11,8 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h"

+#include "GRHydro_Macros.h" + /*@@ @routine GRHydro_Minima_Setup @date Mon Feb 25 11:25:27 2002 @@ -74,17 +76,19 @@ DECLARE_CCTK_PARAMETERS

CCTK_INT i,j,k + CCTK_REAL dummy1 character(len=100) warnline

do i=1,cctk_lsh(1) do j=1,cctk_lsh(2) do k=1,cctk_lsh(3)

- if (rho(i,j,k) < GRHydro_rho_min) then + SET_ATMO_MIN(dummy1, GRHydro_rho_min, r(i,j,k)) + if (rho(i,j,k) < dummy1) then call CCTK_WARN(2,"rho<GRHydro_rho_min!!!") write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel call CCTK_WARN(2,warnline) - write(warnline,'(a25,g15.6)') 'GRHydro_rho_min: ',GRHydro_rho_min + write(warnline,'(a25,g15.6)') 'GRHydro_rho_min: ', dummy1 call CCTK_WARN(2,warnline) write(warnline,'(a25,g15.6)') 'rho: ',rho(i,j,k) call CCTK_WARN(2,warnline)

File [modified]: GRHydro_Reconstruct.F90 Delta lines: +6 -4 =================================================================== --- trunk/src/GRHydro_Reconstruct.F90 2013-01-14 14:23:24 UTC (rev 452) +++ trunk/src/GRHydro_Reconstruct.F90 2013-01-14 14:23:27 UTC (rev 453) @@ -12,6 +12,7 @@ #include "cctk_Parameters.h" #include "cctk_Functions.h"

+#include "GRHydro_Macros.h" #include "SpaceMask.h"

/*@@ @@ -45,7 +46,7 @@ DECLARE_CCTK_FUNCTIONS

CCTK_INT :: i,j,k - CCTK_REAL :: local_min_tracer + CCTK_REAL :: local_min_tracer, dummy1, dummy2 CCTK_INT :: type_bits, not_trivial CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w

@@ -117,12 +118,13 @@ call CCTK_WARN(0, "Flux direction not x,y,z") end if

- !$OMP PARALLEL DO PRIVATE(i,j,k, agxx, agxy, agxz, agyy, agyz, agzz, w) + !$OMP PARALLEL DO PRIVATE(i,j,k, agxx, agxy, agxz, agyy, agyz, agzz, w, dummy1, dummy2) 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 - if(rhoplus(i,j,k).lt.GRHydro_rho_min .or. & - rhominus(i,j,k).lt.GRHydro_rho_min) then + SET_ATMO_MIN(dummy2, GRHydro_rho_min, r(i,j,k)) + if(rhoplus(i,j,k).lt.dummy2 .or. & + rhominus(i,j,k).lt.dummy2) then

rhoplus(i,j,k) = rho(i,j,k) rhominus(i,j,k) = rho(i,j,k)

File [modified]: GRHydro_UpdateMask.F90 Delta lines: +18 -14 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2013-01-14 14:23:24 UTC (rev 452) +++ trunk/src/GRHydro_UpdateMask.F90 2013-01-14 14:23:27 UTC (rev 453) @@ -198,14 +198,17 @@ implicit none

DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS

CCTK_INT :: i,j,k - - !$OMP PARALLEL DO PRIVATE(i,j,k) + CCTK_REAL :: dummy1, dummy2 + + !$OMP PARALLEL DO PRIVATE(i,j,k, dummy1,dummy2) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if (rho(i,j,k) .le. GRHydro_rho_min) then + !if (rho(i,j,k) .le. GRHydro_rho_min) then + IF_BELOW_ATMO(rho(i,j,k), GRHydro_rho_min, 0.0, r(i,j,k)) then atmosphere_mask(i,j,k) = 1 atmosphere_mask_real(i,j,k) = 1 end if @@ -247,7 +250,7 @@ DECLARE_CCTK_PARAMETERS

CCTK_INT :: i, j, k - CCTK_REAL :: det + CCTK_REAL :: det, dummy1, dummy2

! save memory when MP is not used @@ -283,14 +286,14 @@

if (verbose.eq.1) call CCTK_INFO("Entering AtmosphereReset.")

-!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) 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 + SET_ATMO_MIN(rho(i,j,k), GRHydro_rho_min, r(i,j,k)) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.0d0 @@ -359,7 +362,7 @@

CCTK_INT :: i, j, k CCTK_REAL :: det - CCTK_REAL :: rho_min + CCTK_REAL :: rho_min, dummy1, dummy2

CCTK_INT :: eos_handle

@@ -447,15 +450,16 @@ rho_min = rho_min * initial_atmosphere_factor endif

-!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr) +!$OMP PARALLEL DO PRIVATE(det,keytemp,i,j,k,anyerr,keyerr, dummy1, dummy2) 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. & + SET_ATMO_MIN(dummy2, rho_min, r(i,j,k)) + if (rho(i,j,k) .le. dummy2 .or. & GRHydro_enable_internal_excision /= 0 .and. & hydro_excision_mask(i,j,k) .ne. 0) then - rho(i,j,k) = rho_min + SET_ATMO_MIN(rho(i,j,k), dummy2, r(i,j,k)) 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))

@@ -496,8 +500,8 @@ endif end if if (timelevels .gt. 1) then - if (rho_p(i,j,k) .le. rho_min) then - rho_p(i,j,k) = rho_min + if (rho_p(i,j,k) .le. dummy2) then + SET_ATMO_MIN(rho_p(i,j,k), dummy2, r(i,j,k)) velx_p(i,j,k) = 0.0d0 vely_p(i,j,k) = 0.0d0 velz_p(i,j,k) = 0.0d0 @@ -540,8 +544,8 @@ endif 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 + if (rho_p_p(i,j,k) .le. dummy2) then + SET_ATMO_MIN(rho_p_p(i,j,k), dummy2, r(i,j,k)) velx_p_p(i,j,k) = 0.0d0 vely_p_p(i,j,k) = 0.0d0 velz_p_p(i,j,k) = 0.0d0

Directory: /trunk/ ==================

File [modified]: param.ccl Delta lines: +20 -0 =================================================================== --- trunk/param.ccl 2013-01-14 14:23:24 UTC (rev 452) +++ trunk/param.ccl 2013-01-14 14:23:27 UTC (rev 453) @@ -455,6 +455,26 @@ 0.0: :: "Zero or larger. A useful value could be 0.0001" } 0.0

+REAL atmo_falloff_radius "The radius for which the atmosphere starts to falloff as (atmo_falloff_radius/r)**atmo_falloff_power" +{ + 0:* :: "Anything positive" +} 50.0 + +REAL atmo_falloff_power "The power at which the atmosphere level falls off as (atmo_falloof_radius/r)**atmo_falloff_power" +{ + 0:* :: "Anything positive" +} 0.0 + +REAL atmo_tolerance_radius "The radius for which the atmosphere tolerance starts to increase as (r/atmo_tolerance_radius)**atmo_tolerance_power" +{ + 0:* :: "Anything positive" +} 50.0 + +REAL atmo_tolerance_power "The power at which the atmosphere tolerance increases as (r/atmo_tolerance_radius)**atmo_tolerance_power" +{ + 0:* :: "Anything positive" +} 0.0 + BOOLEAN wk_atmosphere "Use some of Wolfgang Kastauns atmosphere tricks" { } "no"