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@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" +}

-