aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_PPM.F90
diff options
context:
space:
mode:
Diffstat (limited to 'src/GRHydro_PPM.F90')
-rw-r--r--src/GRHydro_PPM.F90519
1 files changed, 519 insertions, 0 deletions
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90
index 50da41c..a6c46dc 100644
--- a/src/GRHydro_PPM.F90
+++ b/src/GRHydro_PPM.F90
@@ -1026,6 +1026,525 @@ return
end subroutine SimplePPM_1d
+!!!! routine for doing PPM to temperature
+subroutine SimplePPM_temperature_1d(&
+ nx,dx,velx,temperature,press,&
+ tempminus,&
+ tempplus, hydro_excision_mask)
+
+ USE GRHydro_Scalars
+ USE GRHydro_Eigenproblem
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL, parameter :: one = 1.0d0
+
+ CCTK_INT :: nx
+ CCTK_REAL :: dx
+ CCTK_REAL, dimension(nx) :: temperature, velx
+ CCTK_REAL, dimension(nx) :: tempminus
+ CCTK_REAL, dimension(nx) :: tempplus
+ CCTK_REAL, dimension(nx) :: tempminusl
+ CCTK_REAL, dimension(nx) :: tempplusl
+ CCTK_REAL, dimension(nx) :: tempminusr
+ CCTK_REAL, dimension(nx) :: tempplusr
+ CCTK_REAL, dimension(nx) :: atmosphere_mask
+
+ CCTK_INT :: i,s
+ CCTK_REAL, dimension(nx) :: dtemp
+ CCTK_REAL, dimension(nx) :: dmtemp
+ CCTK_REAL, dimension(nx) :: press,dpress,tilde_flatten
+ CCTK_REAL :: dpress2,dvel,w,flatten,eta,etatilde
+
+ CCTK_INT, dimension(nx) :: hydro_excision_mask
+
+
+ CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax
+
+ logical :: cond
+
+
+#define STEEP(x,dx,dmx) \
+ 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.d0 * abs(x(i) - x(i-1)), \
+ 2.d0 * abs(x(i+1) - x(i))) &&\
+ else &&\
+ dmx(i) = 0.d0 &&\
+ end if
+
+
+ if (use_enhanced_ppm .eq. 0) then
+ !! This is the original PPM algorithm by Colella & Woodward 1984.
+
+ !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178
+ !!$ This is the expression for an even grid.
+
+ do i = 2, nx - 1
+ dpress(i) = press(i+1) - press(i-1)
+ ! the denominator is not necessary
+ dtemp(i) = 0.5d0 * (temperature(i+1) - temperature(i-1))
+ end do
+
+ !!$ Steepened slope. See (1.8) of Colella and Woodward, p.178
+
+ do i = 2, nx - 1
+ STEEP(temperature, dtemp, dmtemp)
+ enddo
+
+
+!!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178
+
+ do i = 2, nx-2
+ tempplus(i) = 0.5d0 * (temperature(i) + temperature(i+1)) + &
+ (dmtemp(i) - dmtemp(i+1)) / 6.d0
+ tempminus(i+1) = tempplus(i)
+ enddo
+ else
+ !! This is the modified PPM algorithm by Colella & Sekora 2008 and McCorquodale & Colella 2011.
+ !! This uses a better limiter based on second derivatives that preserves
+ !! accuracy at local extrema. It also uses a higher-order interpolation polynomial.
+
+ !!$ Initial boundary states (sixth order accurate). See (17) of Colella and Sekora 2008, p.7071
+#define APPROX_AT_CELL_INTERFACE_STENCIL4T(a, ah) \
+ ah = 37.0d0/60.0d0*(a(i)+a(i+1)) - 2.0d0/15.0d0*(a(i-1)+a(i+2)) + 1.0d0/60.0d0*(a(i-2)+a(i+3))
+
+ !!$ Initial boundary states (4th order accurate). See (16) of Colella and Sekora 2008, p.7071
+#define APPROX_AT_CELL_INTERFACET(a, ah) \
+ ah = 7.0d0/12.0d0*(a(i)+a(i+1)) - 1.0d0/12.0d0*(a(i-1)+a(i+2))
+
+#define LIMITT(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) then &&\
+ alim = 0.5d0*(a(i)+a(i+1)) - D2aLim &&\
+ else &&\
+ alim = 0.5d0*(a(i)+a(i+1)) &&\
+ end if &&\
+ endif
+
+#define LIMITT_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
+
+
+
+ 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.
+ do i = 2, nx-2
+
+ APPROX_AT_CELL_INTERFACET(temperature, tempplus(i))
+ LIMITT_EPS(temperature, tempplus(i), enhanced_ppm_C2, tempplus(i))
+ tempminus(i+1) = tempplus(i)
+
+ enddo
+
+ else
+
+ !! Same as above but for 4 stencil points using (17) of Colella & Sekora 2008 as
+ !! initial states.
+ do i = 3, nx-3
+
+ APPROX_AT_CELL_INTERFACE_STENCIL4T(temperature, tempplus(i))
+ LIMITT_EPS(temperature, tempplus(i), enhanced_ppm_C2, tempplus(i))
+ tempminus(i+1) = tempplus(i)
+
+ enddo
+
+ !! Finally compute pressure gradient needed for flattening and shock detection
+ do i = 2, nx-1
+ dpress(i) = press(i+1) - press(i-1)
+ end do
+
+ endif
+ endif
+
+!!$ Zone flattening. See appendix of C&W, p. 197-8.
+ do i = 3, nx - 2
+ dpress2 = press(i+2) - press(i-2)
+ dvel = velx(i+1) - velx(i-1)
+ if ( (abs(dpress(i)) > ppm_epsilon * min(press(i-1),press(i+1))) .and. &
+ (dvel < 0.d0) ) then
+ w = 1.d0
+ else
+ w = 0.d0
+ end if
+ if (abs(dpress2) < ppm_small) then
+ tilde_flatten(i) = 1.d0
+ else
+ tilde_flatten(i) = max(0.d0, 1.d0 - w * max(0.d0, ppm_omega2 * &
+ (dpress(i) / dpress2 - ppm_omega1)))
+ end if
+ end do
+
+
+ if (use_enhanced_ppm .eq. 0) then
+ ! In 1984 PPM, flattening is applied before constraining 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)
+ tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i)
+ tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i)
+ end do
+ else !!$ Really implement C&W, page 197; which requires stencil 4.
+ do i = 4, nx - 3
+ s=sign(one, -dpress(i))
+ flatten = max(tilde_flatten(i), tilde_flatten(i+s))
+ tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i)
+ tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i)
+ end do
+ end if
+ endif
+
+!!$ Monotonicity. See (1.10) of C&W.
+#define MONT(xminus,x,xplus) \
+ if (.not.( (xplus(i).eq.x(i)) .and. (x(i).eq.xminus(i)) ) \
+ .and. ((xplus(i)-x(i))*(x(i)-xminus(i)) .le. 0.d0)) then&&\
+ xminus(i) = x(i) &&\
+ xplus(i) = x(i) &&\
+ else if (6.d0 * (xplus(i) - xminus(i)) * (x(i) - 0.5d0 * \
+ (xplus(i) + xminus(i))) > \
+ (xplus(i) - xminus(i))**2) then &&\
+ xminus(i) = 3.d0 * x(i) - 2.d0 * xplus(i) &&\
+ else if (6.d0 * (xplus(i) - xminus(i)) * (x(i) - 0.5d0 * \
+ (xplus(i) + xminus(i))) < \
+ -(xplus(i) - xminus(i))**2) then &&\
+ xplus(i) = 3.d0 * x(i) - 2.d0 * xminus(i) &&\
+ end if &&\
+
+
+!! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34
+!! This requires 4 stencil points.
+#define MONT_WITH_LOCAL_EXTREMUM_STENCIL4(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 (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 &&\
+ 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 &&\
+ 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 MONT_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 &&\
+ 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 &&\
+ 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
+
+
+!! 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 MONT_WITH_LOCAL_EXTREMUM(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 (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 &&\
+ 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 &&\
+ 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 MONT_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 &&\
+ 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) .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
+
+
+ if (use_enhanced_ppm .eq. 1) then
+ !! Constrain parabolic profiles, PPM 2011/2008
+ if (PPM3) then
+ do i = 3, nx - 2
+ MONT_WITH_LOCAL_EXTREMUM_EPS(tempminus,temperature,tempplus, enhanced_ppm_C2)
+ enddo
+ else
+ do i = 4, nx - 3
+ MONT_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(tempminus,temperature,tempplus, enhanced_ppm_C2, enhanced_ppm_C3)
+ enddo
+ 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)
+ tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i)
+ tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i)
+ end do
+ else !!$ Really implement C&W, page 197; which requires stencil 4.
+ do i = 4, nx - 3
+ s=sign(one, -dpress(i))
+ flatten = max(tilde_flatten(i), tilde_flatten(i+s))
+ tempplus(i) = flatten * tempplus(i) + (1.d0 - flatten) * temperature(i)
+ tempminus(i) = flatten * tempminus(i) + (1.d0 - flatten) * temperature(i)
+ end do
+ end if
+
+ endif
+
+ if (use_enhanced_ppm .eq. 0) then
+ !! Constrain parabolic profiles, PPM 1984
+ do i = GRHydro_stencil, nx - GRHydro_stencil + 1
+ ! original Colella&Woodward monotonicity preservation
+ MONT(tempminus,temperature,tempplus)
+ enddo
+ endif
+
+ !!$ excision
+ do i = 1, nx
+ if (GRHydro_enable_internal_excision /= 0 .and. &
+ (hydro_excision_mask(i) .ne. 0)) then
+ ! do nothing
+ else
+ !!$ Do not optimize cond away by combining the 'if's. Fortran does not
+ !!$ have to follow the order of sub-expressions given here and might
+ !!$ access outside the array range
+ cond = .false.
+ if (i .gt. 1 .and. GRHydro_enable_internal_excision /= 0) then
+ cond = hydro_excision_mask(i-1) .ne. 0
+ end if
+ if (cond) then
+ tempminus(i) = temperature(i)
+ tempplus(i) = temperature(i)
+ tempminus(i-1) = temperature(i)
+ tempplus(i-1) = temperature(i)
+ else
+ cond = .false.
+ if ((i.gt.2) .and. (i.lt.nx) .and. GRHydro_enable_internal_excision /= 0) then
+ cond = (ppm_mppm .eq. 0) .and. (hydro_excision_mask(i-2) .ne. 0)
+ end if
+ if (cond) then
+ call PPM_TVD(temperature(i-1), temperature(i), temperature(i+1), tempminus(i), tempplus(i))
+ end if
+ end if
+ cond = .false.
+ if (i .lt. nx .and. GRHydro_enable_internal_excision /= 0) then
+ cond = hydro_excision_mask(i+1) .ne. 0
+ end if
+ if (cond) then
+ tempminus(i) = temperature(i)
+ tempplus(i) = temperature(i)
+ tempminus(i-1) = temperature(i)
+ tempplus(i-1) = temperature(i)
+ else
+ cond = .false.
+ if ((i.lt.nx-1) .and. (i.gt.1) .and. GRHydro_enable_internal_excision /= 0) then
+ cond = (ppm_mppm .eq. 0) .and. (hydro_excision_mask(i+2) .ne. 0)
+ end if
+ if (cond) then
+ call PPM_TVD(temperature(i-1), temperature(i), temperature(i+1), tempminus(i), tempplus(i))
+ end if
+ end if
+ end if
+ end do
+
+
+
+return
+
+end subroutine SimplePPM_temperature_1d
+
+
+
subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, &
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-07 13:07:39 +0000