User: reisswig
Date: 2012/02/03 05:09 PM
Modified:
/trunk/src/
GRHydro_Eigenproblem.F90, GRHydro_PPM.F90
Log:
Set cs2 to zero if cs2 < 0 and abs(cs2) < 1e-4
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_Eigenproblem.F90
Delta lines: +8 -1
===================================================================
--- trunk/src/GRHydro_Eigenproblem.F90 2012-01-25 01:20:26 UTC (rev 311)
+++ trunk/src/GRHydro_Eigenproblem.F90 2012-02-03 23:09:46 UTC (rev 312)
@@ -57,6 +57,7 @@
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 @@
(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
File [modified]: GRHydro_PPM.F90
Delta lines: +243 -133
===================================================================
--- trunk/src/GRHydro_PPM.F90 2012-01-25 01:20:26 UTC (rev 311)
+++ trunk/src/GRHydro_PPM.F90 2012-02-03 23:09:46 UTC (rev 312)
@@ -99,6 +99,7 @@
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,8 +109,8 @@
!!$ 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 &&\
dmx(i) = sign(1.d0, dx(i)) * \
@@ -195,6 +196,8 @@
!! 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,14 +205,19 @@
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
!! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE,
@@ -518,6 +526,7 @@
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 @@
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 @@
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
+ !! 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
+ 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
- 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)
+ 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
- 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
+ else !!$ Really implement C&W, page 197; which requires stencil 4.
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)
+ 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
- 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.
- 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
+ end if
- endif
endif
if (use_enhanced_ppm .eq. 0) then
@@ -1122,7 +1193,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) 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 @@
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,28 +1306,47 @@
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
!! 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!
@@ -1278,27 +1369,46 @@
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
