diff options
| -rw-r--r-- | param.ccl | 5 | ||||
| -rw-r--r-- | src/GRHydro_PPM.F90 | 227 | ||||
| -rw-r--r-- | src/GRHydro_PPMReconstruct_drv.F90 | 13 |
3 files changed, 138 insertions, 107 deletions
@@ -208,11 +208,6 @@ real enhanced_ppm_C3 "Parameter for enhancecd ppm limiter. Default from McCorquo 0:* :: "must be greater than 0." } 0.1 -real GRHydro_ppm_atmo_tolerance "A relative value for rho below which we switch to standard ppm. This may be necessary for very sharp transitions in the density to atmosphere values." STEERABLE=RECOVER -{ - 0:* :: "anything positive" -} 1e4 - int eno_order "The order of accuracy of the ENO reconstruction" { diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index 6e9c37d..b58f62a 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -57,7 +57,7 @@ subroutine SimplePPM_1d(handle,poly,& velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,& velzplus,epsplus,trivial_rp, hydro_excision_mask,& gxx, gxy, gxz, gyy, gyz, gzz, beta, alp, w_lorentz, & - dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw, rho_min) + dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw) USE GRHydro_Scalars USE GRHydro_Eigenproblem @@ -101,23 +101,21 @@ subroutine SimplePPM_1d(handle,poly,& CCTK_REAL :: dupw, dloc, delta CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz CCTK_REAL, dimension(nx) :: xwind, l_ev_l, l_ev_r - !CCTK_INT, dimension(nx) :: use_std_ppm - CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, threshold, rho_min, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax + CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax logical :: cond !!$ Initially, all the Riemann problems will be trivial - trivial_rp = .true. -!use_std_ppm = 1 + #define STEEP(x,dx,dmx) \ - if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0 ) then &&\ + if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0.d0 ) then &&\ dmx(i) = sign(one, dx(i)) * \ - min(abs(dx(i)), 2 * abs(x(i) - x(i-1)), \ - 2 * abs(x(i+1) - x(i))) &&\ + min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \ + 2.d0 * abs(x(i+1) - x(i))) &&\ else &&\ dmx(i) = 0.d0 &&\ end if @@ -195,11 +193,6 @@ trivial_rp = .true. #define APPROX_AT_CELL_INTERFACE(a, ah) \ ah = 7.0d0/12.0d0*(a(i)+a(i+1)) - 1.0d0/12.0d0*(a(i-1)+a(i+2)) - !! A threshold value of rho - threshold = rho_min * GRHydro_ppm_atmo_tolerance - -! .and. rho(i-2) .gt. threshold .and. rho(i+3) .gt. threshold - #define LIMIT(a,ah,C, alim) \ if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ alim = ah &&\ @@ -208,24 +201,34 @@ trivial_rp = .true. D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ else &&\ alim = 0.5d0*(a(i)+a(i+1)) &&\ end if &&\ endif -! .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1))) +#define LIMIT_EPS(a,ah,C, alim) \ + if ((min(a(i),a(i+1)) .le. ah) .and. (ah .le. max(a(i),a(i+1)))) then &&\ + alim = ah &&\ + else &&\ + D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ + D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ + D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ + D2aLim = sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ + alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\ + else &&\ + alim = 0.5d0*(a(i)+a(i+1)) &&\ + end if &&\ + endif -! .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1))) -!sign(1.0d0, D2a)*min(abs(D2aLim), abs(0.5*(a(i)+a(i+1)))) if (PPM3) then !! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, !! then checking for (13) of Colella & Sekora 2008 and applying !! (18) and (19) if (13) is not satisfied. This is done with LIMIT. - !! There, we also switch to lower order if we are near atmosphere values. do i = 2, nx-2 APPROX_AT_CELL_INTERFACE(rho, rhoplus(i)) LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i)) @@ -246,7 +249,7 @@ trivial_rp = .true. if (poly .eq. 0) then do i = 2, nx-2 APPROX_AT_CELL_INTERFACE(eps, epsplus(i)) - LIMIT(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) + LIMIT_EPS(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) epsminus(i+1) = epsplus(i) end do end if @@ -559,6 +562,60 @@ trivial_rp = .true. D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ + else &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ + endif + +!! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. +!! This contains an additional limiter in case the correction becomes larger than +!! the corrected value. This is to avoid negative values for epsilon. +!! This requires 4 stencil points. +#define MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(aminus, a, aplus, C, C3) \ + daplus = aplus(i)-a(i) &&\ + daminus = a(i)-aminus(i) &&\ + if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aLL = a(i-3) - 2.0d0*a(i-2) + a(i-1) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + D2aRR = a(i+1) - 2.0d0*a(i+2) + a(i+3) &&\ + D3a = D2aR - D2aC &&\ + D3aL = D2aC - D2aL &&\ + D3aR = D2aRR - D2aR &&\ + D3aLL = D2aL - D2aLL &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ + D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ + if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ if (daplus*daminus .lt. 0) then &&\ aplus(i) = a(i) + daplus * rhi &&\ @@ -599,35 +656,8 @@ trivial_rp = .true. trivial_rp(i-1) = .false. &&\ trivial_rp(i) = .false. &&\ endif - -! We use (20) of Colella & Sekora 2008 to check wether we are at a local maximum. -! If so, we make use of (22) and (23) to preserve accuracy there. Otherwise we use -! enhanced monotocinity preservation (26) -! #define MON_WITH_LOCAL_EXTREMUM(aminus, a, aplus, C) \ -! if ((aplus(i)-a(i))*(a(i)-aminus(i)) .le. 0 .or. (a(i-1)-a(i))*(a(i)-a(i+1)) .le. 0) then &&\ -! D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ -! D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ -! D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ -! D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ -! if (D2a .ne. 0 .and. sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ -! aplus(i) = a(i) + (aplus(i) -a(i)) * sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a &&\ -! aminus(i) = a(i) + (aminus(i)-a(i)) * sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a &&\ -! else &&\ -! aplus(i) = a(i) &&\ -! aminus(i) = a(i) &&\ -! end if &&\ -! trivial_rp(i-1) = .false. &&\ -! trivial_rp(i) = .false. &&\ -! else &&\ -! if (abs(aplus(i)-a(i)) .ge. 2.0d0*abs(aminus(i)-a(i))) then &&\ -! aplus(i) = a(i) - 2.0d0*(aminus(i) - a(i)) &&\ -! else if (abs(aminus(i)-a(i)) .ge. 2.0d0*abs(aplus(i)-a(i))) then &&\ -! aminus(i) = a(i) - 2.0d0*(aplus(i) - a(i)) &&\ -! end if &&\ -! trivial_rp(i-1) = .false. &&\ -! trivial_rp(i) = .false. &&\ -! endif + !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. !! This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points! @@ -650,6 +680,51 @@ trivial_rp = .true. rhi = D2aLim / D2a &&\ endif &&\ if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ + if (daplus*daminus .lt. 0) then &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*(1.0d0-rhi)*daplus - rhi*daminus &&\ + else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ + endif &&\ + endif &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ + else &&\ + if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ + aplus(i) = a(i) + 2.0d0*daminus &&\ + else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ + aminus(i) = a(i) - 2.0d0*daplus &&\ + end if &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ + endif + + +!! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. +!! This does not use the check for deviations from a cubic. Thus, it gets away with only 3 stencil points! +!! This contains an additional limiter in case the correction becomes larger than +!! the corrected value. This is to avoid negative values for epsilon. +#define MON_WITH_LOCAL_EXTREMUM_EPS(aminus, a, aplus, C) \ + daplus = aplus(i)-a(i) &&\ + daminus = a(i)-aminus(i) &&\ + if (daplus*daminus .le. 0 .or. (a(i-2)-a(i))*(a(i)-a(i+2)) .le. 0) then &&\ + D2a = - (12.0d0*a(i) - 6.0d0*(aminus(i)+aplus(i))) && \ + D2aC = a(i-1) - 2.0d0*a(i) + a(i+1) &&\ + D2aL = a(i-2) - 2.0d0*a(i-1) + a(i) &&\ + D2aR = a(i) - 2.0d0*a(i+1) + a(i+2) &&\ + if (sign(one, D2a) .eq. sign(one, D2aC) .and. sign(one, D2a) .eq. sign(one, D2aL) .and. sign(one, D2a) .eq. sign(one, D2aR)) then &&\ + D2aLim = sign(one, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\ + else &&\ + D2aLim = 0 &&\ + end if &&\ + if (abs(D2a) .le. 1.d-12*max(abs(a(i-2)), abs(a(i-1)), abs(a(i)), abs(a(i+1)), abs(a(i+2)))) then &&\ + rhi = 0 &&\ + else &&\ + rhi = D2aLim / D2a &&\ + endif &&\ + if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ if (daplus*daminus .lt. 0) then &&\ aplus(i) = a(i) + daplus * rhi &&\ @@ -689,8 +764,8 @@ trivial_rp = .true. trivial_rp(i-1) = .false. &&\ trivial_rp(i) = .false. &&\ endif - - + + if (use_enhanced_ppm .eq. 1) then !! Constrain parabolic profiles, PPM 2011/2008 @@ -703,7 +778,7 @@ trivial_rp = .true. end do if (poly .eq. 0) then do i = 3, nx - 2 - MON_WITH_LOCAL_EXTREMUM(epsminus,eps,epsplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM_EPS(epsminus,eps,epsplus, enhanced_ppm_C2) end do end if @@ -716,7 +791,7 @@ trivial_rp = .true. end do if (poly .eq. 0) then do i = 4, nx - 3 - MON_WITH_LOCAL_EXTREMUM_STENCIL4(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3) end do endif end if @@ -744,7 +819,7 @@ trivial_rp = .true. end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(1.d0, -dpress(i)) + s=sign(one, -dpress(i)) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) if (abs(1.d0 - flatten) > 0.d0) then trivial_rp(i-1) = .false. @@ -1037,7 +1112,7 @@ subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, & else !!$ Really implement C&W, page 197; which requires stencil 4. do itracer=1,number_of_tracers do i = 4, nx - 3 - s=sign(1.d0, -dpress(i)) + s=sign(one, -dpress(i)) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) tracerplus(i,itracer) = flatten * tracerplus(i,itracer) + & (1.d0 - flatten) * tracer(i,itracer) @@ -1130,7 +1205,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & CCTK_REAL :: dpress2,w,flatten,dvel CCTK_REAL :: eta, etatilde - CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, threshold, rho_min, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax + CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax do i = 2, nx - 1 @@ -1199,7 +1274,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & D2a = 3.0d0 * (a(i) - 2.0d0*ah + a(i+1)) &&\ D2aL = (a(i-1) - 2.0d0*a(i) + a(i+1)) &&\ D2aR = (a(i) - 2.0d0*a(i+1) + a(i+2)) &&\ - if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0 .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))) then &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) then &&\ alim = 0.5d0*(a(i)+a(i+1)) - sign(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ else &&\ alim = 0.5d0*(a(i)+a(i+1)) &&\ @@ -1270,7 +1345,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(1.d0, -dpress(i)) + s=sign(one, -dpress(i)) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) Y_e_plus(i) = flatten * Y_e_plus(i) + & (1.d0 - flatten) * Y_e(i) @@ -1312,7 +1387,6 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & D3aMin = min(D3aLL, D3aL, D3a, D3aR) &&\ D3aMax = max(D3aLL, D3aL, D3a, D3aR) &&\ if (C3 * max(abs(D3aMin), abs(D3aMax)) .le. D3aMax-D3aMin) then &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ if (daplus*daminus .lt. 0) then &&\ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ @@ -1321,32 +1395,14 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ endif &&\ - else &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) &&\ - aminus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - endif &&\ - endif &&\ endif &&\ endif &&\ else &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ aplus(i) = a(i) + 2.0d0*daminus &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ aminus(i) = a(i) - 2.0d0*daplus &&\ end if &&\ - else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - end if &&\ - endif &&\ endif @@ -1375,7 +1431,6 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & rhi = D2aLim / D2a &&\ endif &&\ if (.not. (rhi .ge. 1.0d0 - 1.d-12)) then &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ if (daplus*daminus .lt. 0) then &&\ aplus(i) = a(i) + daplus * rhi &&\ aminus(i) = a(i) - daminus * rhi &&\ @@ -1384,31 +1439,13 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ aplus(i) = a(i) + 2.0d0*(1.0d0-rhi)*daminus + rhi*daplus &&\ endif &&\ - else &&\ - if (daplus*daminus .lt. 0) then &&\ - aplus(i) = a(i) &&\ - aminus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - else if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - endif &&\ - endif &&\ endif &&\ else &&\ - if (abs(daplus) .le. abs(a(i)) .and. abs(daminus) .le. abs(a(i))) then &&\ if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ aplus(i) = a(i) - 2.0d0*(aminus(i) - a(i)) &&\ else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ aminus(i) = a(i) - 2.0d0*(aplus(i) - a(i)) &&\ end if &&\ - else &&\ - if (abs(daplus) .ge. 2.0d0*abs(daminus)) then &&\ - aplus(i) = a(i) &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) &&\ - end if &&\ - endif &&\ endif @@ -1438,7 +1475,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & end do else !!$ Really implement C&W, page 197; which requires stencil 4. do i = 4, nx - 3 - s=sign(1.d0, -dpress(i)) + s=sign(one, -dpress(i)) flatten = max(tilde_flatten(i), tilde_flatten(i+s)) Y_e_plus(i) = flatten * Y_e_plus(i) + & (1.d0 - flatten) * Y_e(i) diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90 index 3b51e2d..02486d0 100644 --- a/src/GRHydro_PPMReconstruct_drv.F90 +++ b/src/GRHydro_PPMReconstruct_drv.F90 @@ -151,7 +151,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) enddo !$OMP END PARALLEL DO - !!$ PPM starts: if (flux_direction == 1) then !$OMP PARALLEL DO PRIVATE(i, j, k) @@ -168,7 +167,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) w_lorentz(:,j,k), & 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & GRHydro_mppm_eigenvalue_x_right, & - GRHydro_mppm_xwind, GRHydro_rho_min) + GRHydro_mppm_xwind) do i = 1, nx if (trivial_rp(i,j,k)) then SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) @@ -198,7 +197,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & - press(:,j,k), GRHydro_rho_min) + press(:,j,k)) end do end do !$OMP END PARALLEL DO @@ -219,7 +218,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) w_lorentz(j,:,k), & 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & GRHydro_mppm_eigenvalue_y_right, & - GRHydro_mppm_xwind, GRHydro_rho_min) + GRHydro_mppm_xwind) do i = 1, ny if (trivial_rp(j,i,k)) then SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) @@ -249,7 +248,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & velx(j,:,k),vely(j,:,k),velz(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & - press(j,:,k), GRHydro_rho_min) + press(j,:,k)) end do end do !$OMP END PARALLEL DO @@ -270,7 +269,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) w_lorentz(j,k,:), & 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & GRHydro_mppm_eigenvalue_z_right, & - GRHydro_mppm_xwind, GRHydro_rho_min) + GRHydro_mppm_xwind) do i = 1, nz if (trivial_rp(j,k,i)) then SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) @@ -301,7 +300,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velx(j,k,:),vely(j,k,:),velz(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & - press(j,k,:), GRHydro_rho_min) + press(j,k,:)) end do end do !$OMP END PARALLEL DO |