[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 266)
roland.haas at physics.gatech.edu
roland.haas at physics.gatech.edu
Thu Sep 15 11:42:48 CDT 2011
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
More information about the Commits
mailing list