[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 491)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Wed Mar 27 20:46:27 CDT 2013
User: rhaas
Date: 2013/03/27 08:46 PM
Modified:
/trunk/
interface.ccl, param.ccl, schedule.ccl
/trunk/src/
GRHydro_HLLE.F90
Log:
GRHydro: implemented H viscosity for HLLE solver to eliminate the carbuncle instability. Tested with shocktube and TOV. Not yet tested for problem where carbuncles occur.
From: Christian Reisswig <reisswig at scriwalker.(none)>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_HLLE.F90
Delta lines: +168 -3
===================================================================
--- trunk/src/GRHydro_HLLE.F90 2013-03-06 00:13:43 UTC (rev 490)
+++ trunk/src/GRHydro_HLLE.F90 2013-03-28 01:46:26 UTC (rev 491)
@@ -31,6 +31,118 @@
@@*/
+! eta across face in x-direction
+#define etaX(i,j,k, v, c) \
+ (0.5d0*(abs(v(i+1,j,k,0)-v(i,j,k,0)) + abs(c(i+1,j,k)-c(i,j,k))))
+
+! eta across face in x-direction
+#define etaY(i,j,k, v, c) \
+ (0.5d0*(abs(v(i,j+1,k,1)-v(i,j,k,1)) + abs(c(i,j+1,k)-c(i,j,k))))
+
+! eta across face in x-direction
+#define etaZ(i,j,k, v, c) \
+ (0.5d0*(abs(v(i,j,k+1,2)-v(i,j,k,2)) + abs(c(i,j,k+1)-c(i,j,k))))
+
+
+
+subroutine H_viscosity(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
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: i, j, k, nx, ny, nz
+ CCTK_REAL dpdrho,dpdeps, cs2
+ character*256 :: warnline
+
+ ! 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
+
+! 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
+
+ ! 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
+ else
+ g11 => gxx
+ g12 => gxy
+ g13 => gxz
+ g22 => gyy
+ g23 => gyz
+ g33 => gzz
+ vup => vel
+ end if
+
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
+ !$OMP PARALLEL DO PRIVATE(i,j,k,&
+ !$OMP anyerr, keyerr, keytemp,&
+ !$OMP warnline, xpress,xeps,xtemp,xye, dpdrho, dpdeps, cs2)
+ do k = 1, nz
+ do j = 1, ny
+ do i = 1, nx
+
+ ! set to zero initially
+ eos_c(i,j,k) = 0.0d0
+
+ !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
+
+!!$ Set required fluid quantities
+ if (evolve_temper.ne.1) then
+
+ call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,xpress,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDEps(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,dpdeps,keyerr,anyerr)
+
+ call EOS_Omni_DPressByDRho(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,&
+ rho(i,j,k),eps(i,j,k),xtemp,xye,dpdrho,keyerr,anyerr)
+
+ cs2 = (dpdrho + xpress(1) * dpdeps / (rho(i,j,k)**2))/ &
+ (1.0d0 + eps(i,j,k) + xpress(1)/rho(i,j,k))
+
+ eos_c(i,j,k) = sqrt(cs2)
+
+ else
+
+ call EOS_Omni_cs2(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),cs2,keyerr,anyerr)
+
+ eos_c(i,j,k) = sqrt(cs2)
+
+ end if
+
+ end do
+ end do
+ end do
+
+end subroutine H_viscosity
+
+
subroutine GRHydro_HLLE(CCTK_ARGUMENTS)
USE GRHydro_Eigenproblem
@@ -42,6 +154,7 @@
TARGET gxx, gxy, gxz, gyy, gyz, gzz
TARGET betaa, betab, betac
TARGET betax, betay, betaz
+ TARGET lvel, vel
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
@@ -51,17 +164,21 @@
CCTK_REAL, dimension(5) :: consp,consm_i,fplus,fminus,lamplus
CCTK_REAL, dimension(5) :: f1,lamminus
CCTK_REAL, dimension(5) :: qdiff
- CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det
+ CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det, etabar
CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, &
uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r
CCTK_INT :: type_bits, trivial, not_trivial
+ ! sign requires its arguments to be of identical KIND
+ CCTK_REAL, parameter :: one = 1d0
+
! 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
-
+ CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup
+
if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
g11 => gaa
g12 => gab
@@ -72,6 +189,7 @@
beta1 => betaa
beta2 => betab
beta3 => betac
+ vup => lvel
else
g11 => gxx
g12 => gxy
@@ -82,6 +200,7 @@
beta1 => betax
beta2 => betay
beta3 => betaz
+ vup => vel
end if
if (flux_direction == 1) then
@@ -111,7 +230,7 @@
!$OMP avg_det,gxxh,gxyh,gxzh,gyyh,gyzh,gzzh, &
!$OMP uxxh,uxyh,uxzh,uyyh,uyzh,uzzh, &
!$OMP usendh, charmin, charmax, charpm, &
- !$OMP m)
+ !$OMP m,etabar)
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
@@ -381,7 +500,53 @@
else
call CCTK_WARN(0, "Flux direction not x,y,z")
end if
+
+!!$ Compute H viscosity if requested
+ if (apply_H_viscosity .ne. 0) then
+
+ if (flux_direction == 1) then
+ etabar = max(etaX(i,j,k, vup, eos_c),&
+ etaY(i,j,k, vup, eos_c),&
+ etaY(i+1,j,k, vup, eos_c),&
+ etaY(i,j-1,k, vup, eos_c),&
+ etaY(i+1,j-1,k, vup, eos_c),&
+ etaZ(i,j,k, vup, eos_c),&
+ etaZ(i+1,j,k, vup, eos_c),&
+ etaZ(i,j,k-1, vup, eos_c),&
+ etaZ(i+1,j,k-1, vup, eos_c))
+ else if (flux_direction == 2) then
+ etabar = max(etaY(i,j,k, vup, eos_c),&
+ etaX(i,j,k, vup, eos_c),&
+ etaX(i,j+1,k, vup, eos_c),&
+ etaX(i-1,j,k, vup, eos_c),&
+ etaX(i-1,j+1,k, vup, eos_c),&
+ etaZ(i,j,k, vup, eos_c),&
+ etaZ(i,j+1,k, vup, eos_c),&
+ etaZ(i,j,k-1, vup, eos_c),&
+ etaZ(i,j+1,k-1, vup, eos_c))
+ else if (flux_direction == 3) then
+ etabar = max(etaZ(i,j,k, vup, eos_c),&
+ etaX(i,j,k, vup, eos_c),&
+ etaX(i,j,k+1, vup, eos_c),&
+ etaX(i-1,j,k, vup, eos_c),&
+ etaX(i-1,j,k+1, vup, eos_c),&
+ etaY(i,j,k, vup, eos_c),&
+ etaY(i,j,k+1, vup, eos_c),&
+ etaY(i,j-1,k, vup, eos_c),&
+ etaY(i,j-1,k+1, vup, eos_c))
+ else
+ call CCTK_WARN(0, "Flux direction not x,y,z")
+ end if
+
+ ! modify eigenvalues of Roe's matrix by computed H viscosity
+ do m = 1,5
+ lamplus(m) = sign(one, lamplus(m))*max(abs(lamplus(m)), etabar)
+ lamminus(m) = sign(one, lamminus(m))*max(abs(lamminus(m)), etabar)
+ end do
+ endif
+
+
!!$ Find minimum and maximum wavespeeds
charmin = min(0.d0, lamplus(1), lamplus(2), lamplus(3), &
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +6 -0
===================================================================
--- trunk/interface.ccl 2013-03-06 00:13:43 UTC (rev 490)
+++ trunk/interface.ccl 2013-03-28 01:46:26 UTC (rev 491)
@@ -747,3 +747,9 @@
eos_dpdeps_temp, eos_dpdrho_temp
} "Temporaries for the conservative to primitive conversion"
+CCTK_REAL H_viscosity_temps TYPE=GF tags='Prolongation="None" checkpoint="no"'
+{
+ eos_c
+} "Temporaries for H viscosity"
+
+
File [modified]: param.ccl
Delta lines: +4 -0
===================================================================
--- trunk/param.ccl 2013-03-06 00:13:43 UTC (rev 490)
+++ trunk/param.ccl 2013-03-28 01:46:26 UTC (rev 491)
@@ -310,6 +310,10 @@
"normal" :: "Normal numerical viscosity"
} "fast"
+boolean apply_H_viscosity "H viscosity is useful to fix the carbuncle instability seen at strong shocks"
+{
+} no
+
keyword bound "Which boundary condition to use - FIXME"
{
"flat" :: "Zero order extrapolation"
File [modified]: schedule.ccl
Delta lines: +17 -1
===================================================================
--- trunk/schedule.ccl 2013-03-06 00:13:43 UTC (rev 490)
+++ trunk/schedule.ccl 2013-03-28 01:46:26 UTC (rev 491)
@@ -264,6 +264,16 @@
STORAGE: conformal_state
+
+###############################################
+### Storage for the H viscosity temporaries ###
+###############################################
+
+if (apply_H_viscosity)
+{
+ STORAGE: H_viscosity_temps
+}
+
###############################
### Register startup banner ###
###############################
@@ -1547,5 +1557,11 @@
}
+if (apply_H_viscosity) {
+ SCHEDULE H_viscosity in GRHydroRHS BEFORE FluxTerms
+ {
+ LANG: FORTRAN
+ OPTION: LOCAL
+ } "Compute local temporaries for H viscosity"
+}
-
More information about the Commits
mailing list