diff options
| -rw-r--r-- | param.ccl | 2 | ||||
| -rw-r--r-- | src/GRHydro_Eigenproblem.F90 | 2 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 15 | ||||
| -rw-r--r-- | src/GRHydro_WENOReconstruct.F90 | 58 |
4 files changed, 55 insertions, 22 deletions
@@ -265,7 +265,7 @@ int WENO_order "The order of accuracy of the WENO reconstruction" real weno_eps "WENO epsilon parameter" { 0:* :: "small and positive" -} 1e-6 +} 1e-26 keyword riemann_solver "Which Riemann solver to use" { diff --git a/src/GRHydro_Eigenproblem.F90 b/src/GRHydro_Eigenproblem.F90 index 03576dd..69a6232 100644 --- a/src/GRHydro_Eigenproblem.F90 +++ b/src/GRHydro_Eigenproblem.F90 @@ -88,7 +88,7 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps, & if(cs2.lt.0.0d0) then !$OMP CRITICAL if (abs(cs2) .gt. 1.0d-4) then - write(warnline,'(a50,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho + write(warnline,'(a60,6g16.7)') 'abs(cs2), rho, dpdrho, press*dpdeps/rho**2, eps, press/rho: ', abs(cs2), rho, dpdrho, press * dpdeps / (rho**2), eps, press/rho call CCTK_WARN(1,warnline) call CCTK_WARN(1,"cs2 < 0! Check speed of sound calculation!") cs2 = 0.0d0 diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index 6a2fc27..6c59629 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -49,6 +49,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) CCTK_REAL :: local_min_tracer, dummy1, dummy2 CCTK_INT :: type_bits, not_trivial CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w + character*256 :: warnline ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -130,7 +131,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) 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 - + !.or. epsplus(i,j,k) .lt. 0.0d0 .or. epsminus(i,j,k) .lt. 0.0d0) then rhoplus(i,j,k) = rho(i,j,k) rhominus(i,j,k) = rho(i,j,k) epsplus(i,j,k) = eps(i,j,k) @@ -198,9 +199,21 @@ subroutine Reconstruction(CCTK_ARGUMENTS) end where end if end if + + ! Check if epsilon became negative for ideal gas EOS. + if (evolve_temper.eq.0) then + if (epsplus(i,j,k) .lt. 0.0d0) then + epsplus(i,j,k) = eps(i,j,k) + endif + if (epsminus(i,j,k) .lt. 0.0d0) then + epsminus(i,j,k) = eps(i,j,k) + endif + endif + ! Riemann problem might not be trivial anymore!!!! SpaceMask_SetStateBitsF90(space_mask, i-xoffset, j-yoffset, k-zoffset, type_bits, not_trivial) SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, not_trivial) + enddo enddo enddo diff --git a/src/GRHydro_WENOReconstruct.F90 b/src/GRHydro_WENOReconstruct.F90 index e7273a0..fb03f66 100644 --- a/src/GRHydro_WENOReconstruct.F90 +++ b/src/GRHydro_WENOReconstruct.F90 @@ -158,14 +158,15 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & CCTK_INT :: order, nx, i, j, k, r CCTK_REAL, dimension(nx) :: v, vplus, vminus - CCTK_REAL, dimension(order, 1-order:nx+order) :: vdiff CCTK_INT, dimension(nx) :: hydro_excision_mask logical, dimension(nx) :: trivial_rp logical, dimension(nx) :: excise logical :: normal_weno - CCTK_REAL :: large = 1.d10, gamma1, gamma2, gamma3, beta1, beta2, beta3, w1, w2, w3, wbar1, wbar2, wbar3 + CCTK_REAL :: large = 1.d10, gamma1, gamma2, gamma3, beta1, beta2, beta3, vnorm, betanorm + CCTK_REAL :: wplus1, wplus2, wplus3, wbarplus1, wbarplus2, wbarplus3 + CCTK_REAL :: wminus1, wminus2, wminus3, wbarminus1, wbarminus2, wbarminus3 vminus = 0.d0 vplus = 0.d0 @@ -173,12 +174,9 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & excise = .false. trivial_rp = .false. - - !!$ Initialize excision do i = 1, nx if (GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i) .ne. 0)) then - vdiff(1, i) = large * (-1.d0*(1.d0+mod(i,10)))**i trivial_rp(i) = .true. excise(i) = .true. if (i > 1) then @@ -208,6 +206,7 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & if (normal_weno) then !!$ Compute smoothness indicators +!!$ This is from Tchekhovskoy et al 2007 (WHAM code paper). beta1 = beta_shu(1,1)*v(i-2)**2 & + beta_shu(1,2)*v(i-2)*v(i-1) & + beta_shu(1,3)*v(i-1)**2 & @@ -230,25 +229,46 @@ subroutine GRHydro_WENOReconstruct1d(order, nx, v, vminus, vplus, trivial_rp, & + beta_shu(3,6)*v(i+2)**2 - wbar1 = 1.0d0/16.d0 / (weno_eps + beta1)**2 - wbar2 = 5.0d0/8.d0 / (weno_eps + beta2)**2 - wbar3 = 5.0d0/16.d0 / (weno_eps + beta3)**2 + vnorm = (v(i-2)**2 + v(i-1)**2 + v(i)**2 + v(i+1)**2 + v(i+2)**2) + + beta1 = beta1 + 100.0d0*weno_eps*(vnorm + 1.0d0) + beta2 = beta2 + 100.0d0*weno_eps*(vnorm + 1.0d0) + beta3 = beta3 + 100.0d0*weno_eps*(vnorm + 1.0d0) + + betanorm = beta1 + beta2 + beta3 + + beta1 = beta1 / betanorm + beta2 = beta2 / betanorm + beta3 = beta3 / betanorm + + wbarplus1 = 1.0d0/16.0d0 / (weno_eps + beta1)**2 + wbarplus2 = 5.0d0/8.0d0 / (weno_eps + beta2)**2 + wbarplus3 = 5.0d0/16.0d0 / (weno_eps + beta3)**2 - w1 = wbar1 / (wbar1 + wbar2 + wbar3) - w2 = wbar2 / (wbar1 + wbar2 + wbar3) - w3 = wbar3 / (wbar1 + wbar2 + wbar3) + wplus1 = wbarplus1 / (wbarplus1 + wbarplus2 + wbarplus3) + wplus2 = wbarplus2 / (wbarplus1 + wbarplus2 + wbarplus3) + wplus3 = wbarplus3 / (wbarplus1 + wbarplus2 + wbarplus3) + + wbarminus1 = 5.0d0/16.0d0 / (weno_eps + beta1)**2 + wbarminus2 = 5.0d0/8.0d0 / (weno_eps + beta2)**2 + wbarminus3 = 1.0d0/16.0d0 / (weno_eps + beta3)**2 + wminus1 = wbarminus1 / (wbarminus1 + wbarminus2 + wbarminus3) + wminus2 = wbarminus2 / (wbarminus1 + wbarminus2 + wbarminus3) + wminus3 = wbarminus3 / (wbarminus1 + wbarminus2 + wbarminus3) + !!$ Calculate the reconstruction do j = 1, 5 - vplus(i) = vplus(i) + w1 * weno_coeffs(1,j)*v(i-3+j) & - + w2 * weno_coeffs(2,j)*v(i-3+j) & - + w3 * weno_coeffs(3,j)*v(i-3+j) - vminus(i) = vminus(i) + w1 * weno_coeffs(3,6-j)*v(i-3+j) & - + w2 * weno_coeffs(2,6-j)*v(i-3+j) & - + w3 * weno_coeffs(1,6-j)*v(i-3+j) - + vplus(i) = vplus(i) + wplus1 * weno_coeffs(1,j)*v(i-3+j) & + + wplus2 * weno_coeffs(2,j)*v(i-3+j) & + + wplus3 * weno_coeffs(3,j)*v(i-3+j) + vminus(i) = vminus(i) + wminus1 * weno_coeffs(3,6-j)*v(i-3+j) & + + wminus2 * weno_coeffs(2,6-j)*v(i-3+j) & + + wminus3 * weno_coeffs(1,6-j)*v(i-3+j) end do - + !vminus(i) = v(i) + !vplus(i) = v(i) + end if end do |