User: rhaas Date: 2011/09/15 11:42 AM

Modified: /trunk/ param.ccl /trunk/src/ GRHydro_ParamCheck.F90, GRHydro_Tmunu.F90

Log: add option to smoothly dampen Tmunu to vacuum outside of a given radius

File Changes:

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

File [modified]: GRHydro_ParamCheck.F90 Delta lines: +4 -0 =================================================================== --- trunk/src/GRHydro_ParamCheck.F90 2011-08-19 23:58:27 UTC (rev 265) +++ trunk/src/GRHydro_ParamCheck.F90 2011-09-15 16:42:47 UTC (rev 266) @@ -125,5 +125,9 @@ end if end if

+ if (Tmunu_damping_radius_min .gt. Tmunu_damping_radius_max) then + call CCTK_WARN(0, "Minimum damping radius is greater than maximum damping radius!") + end if + end subroutine GRHydro_ParamCheck

File [modified]: GRHydro_Tmunu.F90 Delta lines: +31 -11 =================================================================== --- trunk/src/GRHydro_Tmunu.F90 2011-08-19 23:58:27 UTC (rev 265) +++ trunk/src/GRHydro_Tmunu.F90 2011-09-15 16:42:47 UTC (rev 266) @@ -70,18 +70,22 @@ implicit none

DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS

CCTK_REAL velxlow, velylow, velzlow CCTK_REAL betaxlow, betaylow, betazlow, beta2 CCTK_REAL utlow, uxlow, uylow, uzlow CCTK_REAL rhoenthalpy CCTK_REAL ut,ux,uy,uz,bst,bsx,bsy,bsz,bs2 + CCTK_REAL dampfac CCTK_INT i,j,k

+!!$ Damping factor + dampfac = 1.0

!$OMP PARALLEL DO PRIVATE(i,j,k,velxlow, velylow, velzlow,& !$OMP betaxlow, betaylow, betazlow, beta2, utlow, uxlow, uylow, uzlow,& - !$OMP rhoenthalpy, ut,ux,uy,uz,bst,bsx,bsy,bsz,bs2) + !$OMP rhoenthalpy, ut,ux,uy,uz,bst,bsx,bsy,bsz,bs2,dampfac) do k = 1, cctk_lsh(3) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) @@ -111,21 +115,37 @@ uylow = velylow uzlow = velzlow

+ +!!$ Initialize damping factor + dampfac = 1.0 +!!$ Apply tanh blending for Tmunu. + if ((Tmunu_damping_radius_min .gt. 0) .and. (r(i,j,k) .gt. Tmunu_damping_radius_min)) then + ! 0.5 * (1.0 - tanh(4.0*(x-x0)/sigma0)) + if (r(i,j,k) .lt. Tmunu_damping_radius_max) then + dampfac = 0.5 * (1.0 - tanh(4.0*(r(i,j,k)-0.5*(Tmunu_damping_radius_max+Tmunu_damping_radius_min))/(Tmunu_damping_radius_max-Tmunu_damping_radius_min))) + else + dampfac = 0.0 + continue ! no need to add anything to Tmunu at the current point (it's zero anyway!) + endif + else + dampfac = 1.0 + endif + !!$ Calculate Tmunu (the lower components!).

- eTtt(i,j,k) = eTtt(i,j,k) + rhoenthalpy*utlow**2 + press(i,j,k)*(beta2 - alp(i,j,k)**2) + eTtt(i,j,k) = eTtt(i,j,k) + dampfac * (rhoenthalpy*utlow**2 + press(i,j,k)*(beta2 - alp(i,j,k)**2))

- eTtx(i,j,k) = eTtx(i,j,k) + rhoenthalpy*utlow*uxlow + press(i,j,k)*betaxlow - eTty(i,j,k) = eTty(i,j,k) + rhoenthalpy*utlow*uylow + press(i,j,k)*betaylow - eTtz(i,j,k) = eTtz(i,j,k) + rhoenthalpy*utlow*uzlow + press(i,j,k)*betazlow + eTtx(i,j,k) = eTtx(i,j,k) + dampfac * (rhoenthalpy*utlow*uxlow + press(i,j,k)*betaxlow) + eTty(i,j,k) = eTty(i,j,k) + dampfac * (rhoenthalpy*utlow*uylow + press(i,j,k)*betaylow) + eTtz(i,j,k) = eTtz(i,j,k) + dampfac * (rhoenthalpy*utlow*uzlow + press(i,j,k)*betazlow)

- eTxx(i,j,k) = eTxx(i,j,k) + rhoenthalpy*uxlow**2 + press(i,j,k)*gxx(i,j,k) - eTyy(i,j,k) = eTyy(i,j,k) + rhoenthalpy*uylow**2 + press(i,j,k)*gyy(i,j,k) - eTzz(i,j,k) = eTzz(i,j,k) + rhoenthalpy*uzlow**2 + press(i,j,k)*gzz(i,j,k) + eTxx(i,j,k) = eTxx(i,j,k) + dampfac * (rhoenthalpy*uxlow**2 + press(i,j,k)*gxx(i,j,k)) + eTyy(i,j,k) = eTyy(i,j,k) + dampfac * (rhoenthalpy*uylow**2 + press(i,j,k)*gyy(i,j,k)) + eTzz(i,j,k) = eTzz(i,j,k) + dampfac * (rhoenthalpy*uzlow**2 + press(i,j,k)*gzz(i,j,k))

- eTxy(i,j,k) = eTxy(i,j,k) + rhoenthalpy*uxlow*uylow + press(i,j,k)*gxy(i,j,k) - eTxz(i,j,k) = eTxz(i,j,k) + rhoenthalpy*uxlow*uzlow + press(i,j,k)*gxz(i,j,k) - eTyz(i,j,k) = eTyz(i,j,k) + rhoenthalpy*uylow*uzlow + press(i,j,k)*gyz(i,j,k) + eTxy(i,j,k) = eTxy(i,j,k) + dampfac * (rhoenthalpy*uxlow*uylow + press(i,j,k)*gxy(i,j,k)) + eTxz(i,j,k) = eTxz(i,j,k) + dampfac * (rhoenthalpy*uxlow*uzlow + press(i,j,k)*gxz(i,j,k)) + eTyz(i,j,k) = eTyz(i,j,k) + dampfac * (rhoenthalpy*uylow*uzlow + press(i,j,k)*gyz(i,j,k))

end do end do

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

File [modified]: param.ccl Delta lines: +16 -0 =================================================================== --- trunk/param.ccl 2011-08-19 23:58:27 UTC (rev 265) +++ trunk/param.ccl 2011-09-15 16:42:47 UTC (rev 266) @@ -532,6 +532,22 @@ { } "no"

+############### Parameters for Tmunu damping (useful in atmosphere) ##################### + +CCTK_REAL Tmunu_damping_radius_min "damping radius at which we start to damp with a tanh function" +{ + -1 :: "damping switched off" + 0:* :: "damping radius at which we start to damp" +} -1 + + +CCTK_REAL Tmunu_damping_radius_max "damping radius at which Tmunu becomes 0" +{ + -1 :: "damping switched off" + 0:* :: "greater than minimum radius above" +} -1 + + ############### temporary parameters to be removed once connected issues are fixed.

boolean con2prim_oct_hack "Disregard c2p failures in oct/rotsym90 boundary cells with xyz < 0" STEERABLE=ALWAYS