diff options
Diffstat (limited to 'src/GRHydro_PPM.F90')
| -rw-r--r-- | src/GRHydro_PPM.F90 | 46 |
1 files changed, 28 insertions, 18 deletions
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index 2f65bcb..cbe31a8 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -52,7 +52,7 @@ end subroutine PPM_TVD (iand(mask((i)),(type_bits)).eq.(state_bits)) -subroutine SimplePPM_1d(handle,poly,& +subroutine SimplePPM_1d(apply_enhanced_ppm,handle,poly,& nx,dx,rho,velx,vely,velz,eps,press,rhominus,& velxminus,velyminus,velzminus,epsminus,rhoplus,velxplus,velyplus,& velzplus,epsplus,trivial_rp, hydro_excision_mask,& @@ -68,6 +68,8 @@ subroutine SimplePPM_1d(handle,poly,& CCTK_REAL, parameter :: one = 1 + logical :: apply_enhanced_ppm + CCTK_INT :: handle,poly,nx CCTK_REAL :: dx CCTK_REAL, dimension(nx) :: rho,velx,vely,velz,eps @@ -131,7 +133,7 @@ trivial_rp = .true. end if - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) 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 @@ -310,11 +312,11 @@ trivial_rp = .true. !!$So this is just dropped from eq. (3.2) of C&W. !!$We can get around this by just rescaling the constant k0 (ppm_k0 here). - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then ! Only for 1984 PPM scheme! if (ppm_detect .ne. 0) then - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then ! make sure drho, d2rho and dmrho are computed everywhere! do i = 2, nx-1 drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) @@ -470,7 +472,7 @@ trivial_rp = .true. end do - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) 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 @@ -777,7 +779,7 @@ trivial_rp = .true. endif - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then do i = 3, nx - 2 @@ -852,7 +854,7 @@ trivial_rp = .true. endif - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 1984 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 ! original Colella&Woodward monotonicity preservation @@ -1019,7 +1021,7 @@ return end subroutine SimplePPM_1d !!!! routine for doing PPM to temperature -subroutine SimplePPM_temperature_1d(& +subroutine SimplePPM_temperature_1d(apply_enhanced_ppm,& nx,dx,velx,temperature,press,& tempminus,& tempplus, hydro_excision_mask) @@ -1031,6 +1033,8 @@ subroutine SimplePPM_temperature_1d(& DECLARE_CCTK_PARAMETERS + logical :: apply_enhanced_ppm + CCTK_REAL, parameter :: one = 1.0d0 CCTK_INT :: nx @@ -1063,7 +1067,7 @@ subroutine SimplePPM_temperature_1d(& end if - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) 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 @@ -1165,7 +1169,7 @@ subroutine SimplePPM_temperature_1d(& do i = 2, nx-1 dpress(i) = press(i+1) - press(i-1) end do - + endif !!$ Zone flattening. See appendix of C&W, p. 197-8. @@ -1187,7 +1191,7 @@ subroutine SimplePPM_temperature_1d(& end do - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) 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 @@ -1439,7 +1443,7 @@ subroutine SimplePPM_temperature_1d(& endif - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then do i = 3, nx - 2 @@ -1469,7 +1473,7 @@ subroutine SimplePPM_temperature_1d(& endif - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then !! Constrain parabolic profiles, PPM 1984 do i = GRHydro_stencil, nx - GRHydro_stencil + 1 ! original Colella&Woodward monotonicity preservation @@ -1535,7 +1539,8 @@ end subroutine SimplePPM_temperature_1d -subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, & +subroutine SimplePPM_tracer_1d(apply_enhanced_ppm,& + nx,dx,rho,velx,vely,velz, & tracer,tracerminus,tracerplus,press) USE GRHydro_Scalars @@ -1544,6 +1549,8 @@ subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, & DECLARE_CCTK_PARAMETERS + logical :: apply_enhanced_ppm + CCTK_REAL, parameter :: one = 1 CCTK_INT :: nx @@ -1749,7 +1756,8 @@ subroutine SimplePPM_tracer_1d(nx,dx,rho,velx,vely,velz, & end subroutine SimplePPM_tracer_1d -subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & +subroutine SimplePPM_ye_1d(apply_enhanced_ppm,& + nx,dx,rho,velx,vely,velz, & Y_e,Y_e_minus,Y_e_plus,press) USE GRHydro_Scalars @@ -1758,6 +1766,8 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & DECLARE_CCTK_PARAMETERS + logical :: apply_enhanced_ppm + CCTK_REAL, parameter :: one = 1 CCTK_INT :: nx @@ -1779,7 +1789,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & end do - if (use_enhanced_PPM .eq. 0) then + if (.not.apply_enhanced_PPM) then !!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 !!$ This is the expression for an even grid. @@ -1900,7 +1910,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & end if end do - if (use_enhanced_ppm .eq. 0) then + if (.not.apply_enhanced_ppm) then if (PPM3) then do i = 3, nx - 2 flatten = tilde_flatten(i) @@ -2018,7 +2028,7 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & - if (use_enhanced_ppm .eq. 1) then + if (apply_enhanced_ppm) then ! Constrain parabolic profiles if (PPM3) then do i = GRHydro_stencil, nx - GRHydro_stencil + 1 |