Set cs2 to zero if cs2 < 0 and abs(cs2) < 1e-4

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@312 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

2 files changed, 252 insertions, 135 deletions

diff --git a/src/GRHydro_Eigenproblem.F90 b/src/GRHydro_Eigenproblem.F90
index 35313fa..605f372 100644
--- a/src/GRHydro_Eigenproblem.F90
+++ b/src/GRHydro_Eigenproblem.F90
@@ -57,6 +57,7 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps, &
   CCTK_REAL lam1,lam2,lam3,lamm,lamp,lamm_nobeta,lamp_nobeta
   CCTK_INT handle
   CCTK_REAL dpdrho,dpdeps,press
+  character*256 :: warnline
 
 ! begin EOS Omni vars
   integer :: n,keytemp,anyerr,keyerr(1)
@@ -85,7 +86,13 @@ subroutine eigenvalues(handle,rho,velx,vely,velz,eps, &
        (1.0d0 + eps + press/rho)
 
   if(cs2.lt.0.0d0) then
-     call CCTK_WARN(0,"cs2 < 0! Check speed of sound calculation!")
+     if (abs(cs2) .gt. 1.0d-4) then
+        write(warnline,'(a50,6g16.7)') '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(0,"cs2 < 0! Check speed of sound calculation!")
+     else
+        cs2 = 0
+     endif
   endif
 
 
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90
index 99e0371..aedeff7 100644
--- a/src/GRHydro_PPM.F90
+++ b/src/GRHydro_PPM.F90
@@ -99,6 +99,7 @@ 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
 
@@ -108,7 +109,7 @@ subroutine SimplePPM_1d(handle,poly,&
 !!$  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.d0 ) then           &&\
@@ -195,6 +196,8 @@ trivial_rp = .true.
       !! 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   &&\
@@ -202,13 +205,18 @@ trivial_rp = .true.
 	 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. rho(i) .gt. threshold .and. rho(i+1) .gt. threshold .and. rho(i-1) .gt. threshold .and. rho(i+2) .gt. threshold) then                     &&\
-	    alim = 0.5d0*(a(i)+a(i+1)) - sign(1.0d0, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\
+	 D2aLim = sign(1.0d0, 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)))
+
+! .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
       
@@ -518,6 +526,7 @@ trivial_rp = .true.
       trivial_rp(i) = .false.                                   &&\
     end if
 
+
 !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34
 !! This requires 4 stencil points.
 #define MON_WITH_LOCAL_EXTREMUM_STENCIL4(aminus, a, aplus, C, C3) \
@@ -548,31 +557,75 @@ 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   &&\
-	       else   &&\
-		  aplus(i)  = a(i) + daplus * rhi     &&\
-		  aminus(i) = a(i) - daminus * rhi    &&\
+	       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    &&\
+		  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 &&\
+	       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      &&\
+	 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*(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                                                    &&\
+	 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 &&\
 	 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(1.0d0, D2a) .eq. sign(1.0d0, D2aC) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aL) .and. sign(1.0d0, D2a) .eq. sign(1.0d0, D2aR)) then &&\
+! 	       aplus(i)  = a(i) + (aplus(i) -a(i)) * sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) / D2a                   &&\
+! 	       aminus(i) = a(i) + (aminus(i)-a(i)) * sign(1.0d0, 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!
@@ -595,103 +648,121 @@ 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   &&\
-	    else   &&\
-	       aplus(i)  = a(i) + daplus * rhi     &&\
-	       aminus(i) = a(i) - daminus * rhi    &&\
+	    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    &&\
+	       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 &&\
+	    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      &&\
 	 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*(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  &&\
+	 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 &&\
 	 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
-	 do i = 3, nx - 2
-	       MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2)
-	       MON_WITH_LOCAL_EXTREMUM(velxminus,velx,velxplus, enhanced_ppm_C2)
-	       MON_WITH_LOCAL_EXTREMUM(velyminus,vely,velyplus, enhanced_ppm_C2)
-	       MON_WITH_LOCAL_EXTREMUM(velzminus,velz,velzplus, enhanced_ppm_C2)
-	 end do
-	 if (poly .eq. 0) then
-	       do i = 3, nx - 2
-		  MON_WITH_LOCAL_EXTREMUM(epsminus,eps,epsplus, enhanced_ppm_C2)
-	       end do
-	 end if
-	 
-     else
-	 do i = 4, nx - 3
-	       MON_WITH_LOCAL_EXTREMUM_STENCIL4(rhominus,rho,rhoplus, enhanced_ppm_C2, enhanced_ppm_C3)
-	       MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,velx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3)
-	       MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vely,velyplus, enhanced_ppm_C2, enhanced_ppm_C3)
-	       MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,velz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3)
-	 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)
-	       end do
-	 endif
-
-         ! Apply flattening after constraining the parabolic profiles
-	 if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3.
+      !! Constrain parabolic profiles, PPM 2011/2008
+      if (PPM3) then
 	    do i = 3, nx - 2
-	       flatten = tilde_flatten(i)
-	       if (abs(1.d0 - flatten) > 0.d0) then
-		  trivial_rp(i-1) = .false.
-		  trivial_rp(i) = .false.
-	       end if
-	       rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
-	       rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
-	       velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
-	       velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
-	       velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
-	       velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
-	       velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
-	       velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
-	       if (poly .eq. 0) then
-		  epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
-		  epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
-	       end if
+		  MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2)
+		  MON_WITH_LOCAL_EXTREMUM(velxminus,velx,velxplus, enhanced_ppm_C2)
+		  MON_WITH_LOCAL_EXTREMUM(velyminus,vely,velyplus, enhanced_ppm_C2)
+		  MON_WITH_LOCAL_EXTREMUM(velzminus,velz,velzplus, enhanced_ppm_C2)
 	    end do
-	 else  !!$ Really implement C&W, page 197; which requires stencil 4.
+	    if (poly .eq. 0) then
+		  do i = 3, nx - 2
+		     MON_WITH_LOCAL_EXTREMUM(epsminus,eps,epsplus, enhanced_ppm_C2)
+		  end do
+	    end if
+	    
+      else
 	    do i = 4, nx - 3
-	       s=sign(1.d0, -dpress(i))
-	       flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
-	       if (abs(1.d0 - flatten) > 0.d0) then
-		  trivial_rp(i-1) = .false.
-		  trivial_rp(i) = .false.
-	       end if
-	       rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
-	       rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
-	       velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
-	       velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
-	       velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
-	       velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
-	       velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
-	       velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
-	       if (poly .eq. 0) then
-		  epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
-		  epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
-	       end if
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(rhominus,rho,rhoplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velxminus,velx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vely,velyplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		  MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,velz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3)
 	    end do
-	 end if
+	    if (poly .eq. 0) then
+		  do i = 4, nx - 3
+		     MON_WITH_LOCAL_EXTREMUM_STENCIL4(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3)
+		  end do
+	    endif
+      end if 
+      
+      ! Apply flattening after constraining the parabolic profiles
+      if (PPM3) then !!$ Implement C&W, page 197, but with a workaround which allows to use stencil=3.
+	 do i = 3, nx - 2
+	    flatten = tilde_flatten(i)
+	    if (abs(1.d0 - flatten) > 0.d0) then
+	       trivial_rp(i-1) = .false.
+	       trivial_rp(i) = .false.
+	    end if
+	    rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
+	    rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
+	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
+	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
+	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
+	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
+	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
+	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
+	    if (poly .eq. 0) then
+	       epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
+	       epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
+	    end if
+	 end do
+      else  !!$ Really implement C&W, page 197; which requires stencil 4.
+	 do i = 4, nx - 3
+	    s=sign(1.d0, -dpress(i))
+	    flatten = max(tilde_flatten(i), tilde_flatten(i+s))  
+	    if (abs(1.d0 - flatten) > 0.d0) then
+	       trivial_rp(i-1) = .false.
+	       trivial_rp(i) = .false.
+	    end if
+	    rhoplus(i) = flatten * rhoplus(i) + (1.d0 - flatten) * rho(i)
+	    rhominus(i) = flatten * rhominus(i) + (1.d0 - flatten) * rho(i)
+	    velxplus(i) = flatten * velxplus(i) + (1.d0 - flatten) * velx(i)
+	    velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * velx(i)
+	    velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vely(i)
+	    velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vely(i)
+	    velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * velz(i)
+	    velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * velz(i)
+	    if (poly .eq. 0) then
+	       epsplus(i) = flatten * epsplus(i) + (1.d0 - flatten) * eps(i)
+	       epsminus(i) = flatten * epsminus(i) + (1.d0 - flatten) * eps(i)
+	    end if
+	 end do
+      end if
     
-      endif
   endif
 
   if (use_enhanced_ppm .eq. 0) then
@@ -1122,7 +1193,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) then                     &&\
+	 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)) - sign(1.0d0, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\
 	 else                                                                             &&\
 	    alim = 0.5d0*(a(i)+a(i+1))                                                    &&\
@@ -1203,6 +1274,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
       end if
   endif
   
+
 #undef MON_WITH_LOCAL_EXTREMUM_STENCIL4
 !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34
 !! This requires 4 stencil points.
@@ -1234,26 +1306,45 @@ 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 (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   &&\
-	       else   &&\
-		  aplus(i)  = a(i) + daplus * rhi     &&\
-		  aminus(i) = a(i) - daminus * rhi    &&\
+	       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    &&\
+		  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 &&\
+	       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      &&\
+	 endif &&\
       else                                                         &&\
-	 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                                                    &&\
+	 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
+            
+
+
 
 
 #undef MON_WITH_LOCAL_EXTREMUM
@@ -1278,27 +1369,46 @@ 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 (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   &&\
-	    else   &&\
-	       aplus(i)  = a(i) + daplus * rhi     &&\
-	       aminus(i) = a(i) - daminus * rhi    &&\
+	    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    &&\
+	       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 &&\
+	    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) .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  &&\
+	 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
 
 
+
+
+
   if (use_enhanced_ppm .eq. 1) then
       ! Constrain parabolic profiles
       if (PPM3) then