User: reisswig
Date: 2012/05/14 01:44 PM
Modified:
/trunk/
param.ccl
/trunk/src/
GRHydro_PPM.F90, GRHydro_PPMReconstruct_drv.F90
Log:
Cleanup of PPM routine. Remove unused PPM parameter. Only apply stricter "beyond Colella" limiting scheme to epsilon (to avoid that correction in one step leads to negative values).
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_PPM.F90
Delta lines: +132 -95
===================================================================
--- trunk/src/GRHydro_PPM.F90 2012-05-14 15:37:35 UTC (rev 332)
+++ trunk/src/GRHydro_PPM.F90 2012-05-14 18:44:42 UTC (rev 333)
@@ -57,7 +57,7 @@
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 @@
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 @@
#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,6 +201,21 @@
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) then &&\
+ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\
+ else &&\
+ alim = 0.5d0*(a(i)+a(i+1)) &&\
+ end if &&\
+ endif
+
+#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 &&\
@@ -215,17 +223,12 @@
end if &&\
endif
-! .and. abs(D2aLim) .lt. abs(0.5d0*(a(i)+a(i+1)))
-! .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 @@
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 @@
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,36 +656,9 @@
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!
#define MON_WITH_LOCAL_EXTREMUM(aminus, a, aplus, C) \
@@ -650,6 +680,51 @@
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,9 +764,9 @@
trivial_rp(i-1) = .false. &&\
trivial_rp(i) = .false. &&\
endif
-
-
+
+
if (use_enhanced_ppm .eq. 1) then
!! Constrain parabolic profiles, PPM 2011/2008
if (PPM3) then
@@ -703,7 +778,7 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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)
File [modified]: GRHydro_PPMReconstruct_drv.F90
Delta lines: +6 -7
===================================================================
--- trunk/src/GRHydro_PPMReconstruct_drv.F90 2012-05-14 15:37:35 UTC (rev 332)
+++ trunk/src/GRHydro_PPMReconstruct_drv.F90 2012-05-14 18:44:42 UTC (rev 333)
@@ -151,7 +151,6 @@
enddo
!$OMP END PARALLEL DO
-
!!$ PPM starts:
if (flux_direction == 1) then
!$OMP PARALLEL DO PRIVATE(i, j, k)
@@ -168,7 +167,7 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +0 -5
===================================================================
--- trunk/param.ccl 2012-05-14 15:37:35 UTC (rev 332)
+++ trunk/param.ccl 2012-05-14 18:44:42 UTC (rev 333)
@@ -208,12 +208,7 @@
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"
{
1:* :: "Anything above 1, but above 5 is pointless"
