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.F9046
1 files changed, 26 insertions, 20 deletions
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90
index 6bc5880..6e9c37d 100644
--- a/src/GRHydro_PPM.F90
+++ b/src/GRHydro_PPM.F90
@@ -66,6 +66,8 @@ subroutine SimplePPM_1d(handle,poly,&
DECLARE_CCTK_PARAMETERS
+ CCTK_REAL, parameter :: one = 1
+
CCTK_INT :: handle,poly,nx
CCTK_REAL :: dx
CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps
@@ -112,10 +114,10 @@ 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)) * \
- min(abs(dx(i)), 2.d0 * abs(x(i) - x(i-1)), \
- 2.d0 * abs(x(i+1) - x(i))) &&\
+ if ( (x(i+1) - x(i)) * (x(i) - x(i-1)) > 0 ) 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))) &&\
else &&\
dmx(i) = 0.d0 &&\
end if
@@ -205,7 +207,7 @@ 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)) &&\
- D2aLim = sign(1.0d0, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\
+ 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 &&\
@@ -478,7 +480,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.
@@ -543,8 +545,8 @@ trivial_rp = .true.
D3aL = D2aC - D2aL &&\
D3aR = D2aRR - D2aR &&\
D3aLL = D2aL - D2aLL &&\
- if (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 &&\
- D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\
+ 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 &&\
@@ -608,9 +610,9 @@ trivial_rp = .true.
! 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 &&\
+! 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) &&\
@@ -637,8 +639,8 @@ trivial_rp = .true.
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(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 &&\
- D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\
+ 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 &&\
@@ -898,6 +900,8 @@ subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, &
DECLARE_CCTK_PARAMETERS
+ CCTK_REAL, parameter :: one = 1
+
CCTK_INT :: nx
CCTK_REAL :: dx
CCTK_REAL, dimension(nx) :: rho,velx,vely,velz
@@ -934,7 +938,7 @@ subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, &
do i = 2, nx - 1
if( (tracer(i+1,itracer) - tracer(i,itracer)) * &
(tracer(i,itracer) - tracer(i-1,itracer)) > 0.0d0 ) then
- dmtracer(i,itracer) = sign(1.0d0,dtracer(i,itracer)) * &
+ dmtracer(i,itracer) = sign(one,dtracer(i,itracer)) * &
min(abs(dtracer(i,itracer)), 2.0d0 * &
abs(tracer(i,itracer) - tracer(i-1,itracer)), &
2.0d0 * abs(tracer(i+1,itracer) - tracer(i,itracer)))
@@ -1111,6 +1115,8 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
DECLARE_CCTK_PARAMETERS
+ CCTK_REAL, parameter :: one = 1
+
CCTK_INT :: nx
CCTK_REAL :: dx
CCTK_REAL, dimension(nx) :: rho,velx,vely,velz
@@ -1151,7 +1157,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
do i = 2, nx - 1
if( (Y_e(i+1) - Y_e(i)) * &
(Y_e(i) - Y_e(i-1)) > 0.0d0 ) then
- dmY_e(i) = sign(1.0d0,dY_e(i)) * &
+ dmY_e(i) = sign(one,dY_e(i)) * &
min(abs(dY_e(i)), 2.0d0 * &
abs(Y_e(i) - Y_e(i-1)), &
2.0d0 * abs(Y_e(i+1) - Y_e(i)))
@@ -1194,7 +1200,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
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 &&\
- alim = 0.5d0*(a(i)+a(i+1)) - sign(1.0d0, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\
+ 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)) &&\
end if &&\
@@ -1292,8 +1298,8 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
D3aL = D2aC - D2aL &&\
D3aR = D2aRR - D2aR &&\
D3aLL = D2aL - D2aLL &&\
- if (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 &&\
- D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\
+ 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 &&\
@@ -1358,8 +1364,8 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, &
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(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 &&\
- D2aLim = sign(1.0d0, D2a) * min(C*abs(D2aL), C*abs(D2aR), C*abs(D2aC), abs(D2a)) &&\
+ 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 &&\
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-29 23:25:10 +0000