diff options
| -rw-r--r-- | param.ccl | 21 | ||||
| -rw-r--r-- | src/GRHydro_PPM.F90 | 1121 | ||||
| -rw-r--r-- | src/GRHydro_PPMReconstruct_drv.F90 | 12 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 26 |
4 files changed, 844 insertions, 336 deletions
@@ -193,6 +193,27 @@ int ppm_mppm_debug_eigenvalues "write eigenvalues into debug grid variables" 0:1 :: "0 (off, default) or 1 (on)" } 0 + +boolean use_enhanced_ppm "Use the enhanced ppm reconstruction method proposed by Colella & Sekora 2008 and McCorquodale & Colella 2011" STEERABLE=RECOVER +{ +} no + +real enhanced_ppm_C2 "Parameter for enhancecd ppm limiter. Default from McCorquodale & Colella 2011" +{ + *:* :: "must be greater than 1. According to Colella&Sekora 2008, enhanced ppm is insensitive to C in [1.25,5]" +} 1.25 + +real enhanced_ppm_C3 "Parameter for enhancecd ppm limiter. Default from McCorquodale & Colella 2011" +{ + 0:* :: "must be greater than 0." +} 0.1 + +real GRHydro_ppm_atmo_tolerance "A relative value for rho below which we switch to standard ppm. This may be necessary for very sharp transitions in the density to atmosphere values." STEERABLE=RECOVER +{ + 0:* :: "anything positive" +} 1e4 + + int eno_order "The order of accuracy of the ENO reconstruction" { 1:* :: "Anything above 1, but above 5 is pointless" diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index c5754a0..99e0371 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -32,7 +32,7 @@ end subroutine PPM_TVD /*@@ @routine SimplePPM_1d @date Thu Feb 14 19:08:52 2002 - @author Ian Hawke, Toni Font + @author Ian Hawke, Toni Font, Christian Reisswig @desc The simple PPM reconstruction routine that applies along each one dimensional slice. @@ -43,16 +43,21 @@ end subroutine PPM_TVD @history Written in frustration when IH couldn''t get Toni''s original code to work. + Later extended to enhanced PPM scheme by CR. @endhistory @@*/ +#define SpaceMask_CheckStateBitsF90_1D(mask,i,type_bits,state_bits) \ + (iand(mask((i)),(type_bits)).eq.(state_bits)) + + subroutine SimplePPM_1d(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,& gxx, gxy, gxz, gyy, gyz, gzz, beta, alp, w_lorentz, & - dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw) + dir, ni, nj, nrx, nry, nrz, ev_l, ev_r, xw, rho_min) USE GRHydro_Scalars USE GRHydro_Eigenproblem @@ -72,6 +77,7 @@ subroutine SimplePPM_1d(handle,poly,& CCTK_REAL, dimension(nx) :: rhominusr,velxminusr,velyminusr,velzminusr CCTK_REAL, dimension(nx) :: epsminusr CCTK_REAL, dimension(nx) :: rhoplusr,velxplusr,velyplusr,velzplusr,epsplusr + CCTK_REAL, dimension(nx) :: atmosphere_mask CCTK_INT :: i,s CCTK_REAL, dimension(nx) :: drho,dvelx,dvely,dvelz,deps @@ -94,6 +100,8 @@ subroutine SimplePPM_1d(handle,poly,& CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz CCTK_REAL, dimension(nx) :: xwind, l_ev_l, l_ev_r + CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, threshold, rho_min, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax + logical :: cond @@ -101,70 +109,175 @@ subroutine SimplePPM_1d(handle,poly,& trivial_rp = .true. -!!$ 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 - drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) - dvelx(i) = 0.5d0 * (velx(i+1) - velx(i-1)) - dvely(i) = 0.5d0 * (vely(i+1) - vely(i-1)) - dvelz(i) = 0.5d0 * (velz(i+1) - velz(i-1)) - dpress(i) = press(i+1) - press(i-1) - d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx - ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 - ! the denominator is not necessary - end do - if (poly .eq. 0) then - do i = 2, nx - 1 - deps(i) = 0.5d0 * (eps(i+1) - eps(i-1)) - end do - end if +#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))) &&\ + 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 + drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) + dvelx(i) = 0.5d0 * (velx(i+1) - velx(i-1)) + dvely(i) = 0.5d0 * (vely(i+1) - vely(i-1)) + dvelz(i) = 0.5d0 * (velz(i+1) - velz(i-1)) + dpress(i) = press(i+1) - press(i-1) + d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1))! / 6.d0 / dx / dx + ! since we use d2rho only for the condition d2rho(i+1)*d2rhoi(i-1)<0 + ! the denominator is not necessary + end do + if (poly .eq. 0) then + do i = 2, nx - 1 + deps(i) = 0.5d0 * (eps(i+1) - eps(i-1)) + end do + end if -!!$ Steepened slope. See (1.8) of Colella and Woodward, p.178 + !!$ Steepened slope. See (1.8) of Colella and Woodward, p.178 - do i = 2, nx - 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))) &&\ - else &&\ - dmx(i) = 0.d0 &&\ - end if - STEEP(rho, drho, dmrho) - STEEP(velx, dvelx, dmvelx) - STEEP(vely, dvely, dmvely) - STEEP(velz, dvelz, dmvelz) - end do - if (poly .eq. 0) then - do i = 2, nx - 1 - STEEP(eps, deps, dmeps) - end do - end if + do i = 2, nx - 1 + STEEP(rho, drho, dmrho) + STEEP(velx, dvelx, dmvelx) + STEEP(vely, dvely, dmvely) + STEEP(velz, dvelz, dmvelz) + end do + if (poly .eq. 0) then + do i = 2, nx - 1 + STEEP(eps, deps, dmeps) + end do + end if -!!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178 + !!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178 + + do i = 2, nx-2 + rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + & + (dmrho(i) - dmrho(i+1)) / 6.d0 + rhominus(i+1) = rhoplus(i) + velxplus(i) = 0.5d0 * (velx(i) + velx(i+1)) + & + (dmvelx(i) - dmvelx(i+1)) / 6.d0 + velxminus(i+1) = velxplus(i) + velyplus(i) = 0.5d0 * (vely(i) + vely(i+1)) + & + (dmvely(i) - dmvely(i+1)) / 6.d0 + velyminus(i+1) = velyplus(i) + velzplus(i) = 0.5d0 * (velz(i) + velz(i+1)) + & + (dmvelz(i) - dmvelz(i+1)) / 6.d0 + velzminus(i+1) = velzplus(i) + end do + if (poly .eq. 0) then + do i = 2, nx-2 + epsplus(i) = 0.5d0 * (eps(i) + eps(i+1)) + & + (dmeps(i) - dmeps(i+1)) / 6.d0 + epsminus(i+1) = epsplus(i) + end do + end if + 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_STENCIL4(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_INTERFACE(a, ah) \ + ah = 7.0d0/12.0d0*(a(i)+a(i+1)) - 1.0d0/12.0d0*(a(i-1)+a(i+2)) + + !! A threshold value of rho + threshold = rho_min * GRHydro_ppm_atmo_tolerance + +#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 &&\ + 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)) &&\ + 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 &&\ + 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. + !! There, we also switch to lower order if we are near atmosphere values. + do i = 2, nx-2 + APPROX_AT_CELL_INTERFACE(rho, rhoplus(i)) + LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i)) + rhominus(i+1) = rhoplus(i) + + APPROX_AT_CELL_INTERFACE(velxplus, velxplus(i)) + LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i)) + velxminus(i+1) = velxplus(i) + + APPROX_AT_CELL_INTERFACE(velyplus, velyplus(i)) + LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i)) + velyminus(i+1) = velyplus(i) + + APPROX_AT_CELL_INTERFACE(velzplus, velzplus(i)) + LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i)) + velzminus(i+1) = velzplus(i) + end do + if (poly .eq. 0) then + do i = 2, nx-2 + APPROX_AT_CELL_INTERFACE(eps, epsplus(i)) + LIMIT(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) + epsminus(i+1) = epsplus(i) + end do + end if + + else - do i = 2, nx-2 - rhoplus(i) = 0.5d0 * (rho(i) + rho(i+1)) + & - (dmrho(i) - dmrho(i+1)) / 6.d0 - rhominus(i+1) = rhoplus(i) - velxplus(i) = 0.5d0 * (velx(i) + velx(i+1)) + & - (dmvelx(i) - dmvelx(i+1)) / 6.d0 - velxminus(i+1) = velxplus(i) - velyplus(i) = 0.5d0 * (vely(i) + vely(i+1)) + & - (dmvely(i) - dmvely(i+1)) / 6.d0 - velyminus(i+1) = velyplus(i) - velzplus(i) = 0.5d0 * (velz(i) + velz(i+1)) + & - (dmvelz(i) - dmvelz(i+1)) / 6.d0 - velzminus(i+1) = velzplus(i) - end do - if (poly .eq. 0) then - do i = 2, nx-2 - epsplus(i) = 0.5d0 * (eps(i) + eps(i+1)) + & - (dmeps(i) - dmeps(i+1)) / 6.d0 - epsminus(i+1) = epsplus(i) + !! 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_STENCIL4(rho, rhoplus(i)) + LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i)) + rhominus(i+1) = rhoplus(i) + + APPROX_AT_CELL_INTERFACE_STENCIL4(velxplus, velxplus(i)) + LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i)) + velxminus(i+1) = velxplus(i) + + APPROX_AT_CELL_INTERFACE_STENCIL4(velyplus, velyplus(i)) + LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i)) + velyminus(i+1) = velyplus(i) + + APPROX_AT_CELL_INTERFACE_STENCIL4(velzplus, velzplus(i)) + LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i)) + velzminus(i+1) = velzplus(i) + end do + if (poly .eq. 0) then + do i = 3, nx-3 + APPROX_AT_CELL_INTERFACE_STENCIL4(eps, epsplus(i)) + LIMIT(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) + epsminus(i+1) = epsplus(i) + end do + end if + endif + + !! 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 - end if + + endif + !!$Discontinuity steepening. See (1.14-17) of C&W. !!$This is the detect routine which mat be activated with the ppm_detect parameter @@ -173,135 +286,148 @@ 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 (ppm_detect .ne. 0) then - - do i = 3, nx - 2 - if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - & - ppm_epsilon_shock * min(abs(rho(i+1)), abs(rho(i-1))) > 0.d0) ) then - etatilde = (rho(i-2) - rho(i+2) + 4.d0 * drho(i)) / (drho(i) * 12.d0) - else - etatilde = 0.d0 - end if - eta = max(0.d0, min(1.d0, ppm_eta1 * (etatilde - ppm_eta2))) - if (ppm_k0 * abs(drho(i)) * min(press(i-1),press(i+1)) < & - abs(dpress(i)) * min(rho(i-1), rho(i+1))) then - eta = 0.d0 - end if - if (eta > 0.d0) then - trivial_rp(i-1) = .false. - trivial_rp(i) = .false. - end if - rhominus(i) = rhominus(i) * (1.d0 - eta) + & - (rho(i-1) + 0.5d0 * dmrho(i-1)) * eta - rhoplus(i) = rhoplus(i) * (1.d0 - eta) + & - (rho(i+1) - 0.5d0 * dmrho(i+1)) * eta - end do - - end if - + if (use_enhanced_ppm .eq. 0) then + ! Only for 1984 PPM scheme! + if (ppm_detect .ne. 0) then + + if (use_enhanced_ppm .eq. 1) 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)) + d2rho(i) = (rho(i+1) - 2.d0 * rho(i) + rho(i-1)) + end do + do i = 2, nx-1 + STEEP(rho, drho, dmrho) + end do + endif + + do i = 3, nx - 2 + if ( (d2rho(i+1)*d2rho(i-1) < 0.d0).and.(abs(rho(i+1)-rho(i-1)) - & + ppm_epsilon_shock * min(abs(rho(i+1)), abs(rho(i-1))) > 0.d0) ) then + etatilde = (rho(i-2) - rho(i+2) + 4.d0 * drho(i)) / (drho(i) * 12.d0) + else + etatilde = 0.d0 + end if + eta = max(0.d0, min(1.d0, ppm_eta1 * (etatilde - ppm_eta2))) + if (ppm_k0 * abs(drho(i)) * min(press(i-1),press(i+1)) < & + abs(dpress(i)) * min(rho(i-1), rho(i+1))) then + eta = 0.d0 + end if + if (eta > 0.d0) then + trivial_rp(i-1) = .false. + trivial_rp(i) = .false. + end if + rhominus(i) = rhominus(i) * (1.d0 - eta) + & + (rho(i-1) + 0.5d0 * dmrho(i-1)) * eta + rhoplus(i) = rhoplus(i) * (1.d0 - eta) + & + (rho(i+1) - 0.5d0 * dmrho(i+1)) * eta + end do + + end if + !!$ mppm #define D_UPW(x) (0.5d0 * (x(i) + x(i+1))) #define LEFT1(x) (13.d0*x(i+1)-5.d0*x(i+2)+x(i+3)+3.d0*x(i ))/12.d0 #define RIGHT1(x) (13.d0*x(i )-5.d0*x(i-1)+x(i-2)+3.d0*x(i+1))/12.d0 - if (ppm_mppm .gt. 0) then - l_ev_l=0.d0 - l_ev_r=0.d0 - xwind=0.d0 - do i=3, nx - 3 - agxx = 0.5d0*( gxx(i) + gxx(i+1) ) - agxy = 0.5d0*( gxy(i) + gxy(i+1) ) - agxz = 0.5d0*( gxz(i) + gxz(i+1) ) - agyy = 0.5d0*( gyy(i) + gyy(i+1) ) - agyz = 0.5d0*( gyz(i) + gyz(i+1) ) - agzz = 0.5d0*( gzz(i) + gzz(i+1) ) - det = SPATIAL_DETERMINANT(agxx, agxy, agxz, \ - agyy, agyz, agzz) - call UpperMetric (uxx, uxy, uxz, uyy, uyz, uzz, & - det, agxx, agxy, agxz, agyy, agyz, agzz) - call eigenvalues(handle,& - D_UPW(rho), D_UPW(velx), D_UPW(vely), D_UPW(velz), & - D_UPW(eps), D_UPW(w_lorentz), lam, & - agxx, agxy, agxz, agyy, agyz, agzz, & - uxx, D_UPW(alp), D_UPW(beta)) - l_ev_l(i)=lam(1) - l_ev_r(i)=lam(5) - xwind(i) = (lam(1) + lam(5)) / (abs(lam(1)) + abs(lam(5))) - xwind(i) = min(1.d0, max(-1.d0, xwind(i))) + if (ppm_mppm .gt. 0) then + l_ev_l=0.d0 + l_ev_r=0.d0 + xwind=0.d0 + do i=3, nx - 3 + agxx = 0.5d0*( gxx(i) + gxx(i+1) ) + agxy = 0.5d0*( gxy(i) + gxy(i+1) ) + agxz = 0.5d0*( gxz(i) + gxz(i+1) ) + agyy = 0.5d0*( gyy(i) + gyy(i+1) ) + agyz = 0.5d0*( gyz(i) + gyz(i+1) ) + agzz = 0.5d0*( gzz(i) + gzz(i+1) ) + det = SPATIAL_DETERMINANT(agxx, agxy, agxz, \ + agyy, agyz, agzz) + call UpperMetric (uxx, uxy, uxz, uyy, uyz, uzz, & + det, agxx, agxy, agxz, agyy, agyz, agzz) + call eigenvalues(handle,& + D_UPW(rho), D_UPW(velx), D_UPW(vely), D_UPW(velz), & + D_UPW(eps), D_UPW(w_lorentz), lam, & + agxx, agxy, agxz, agyy, agyz, agzz, & + uxx, D_UPW(alp), D_UPW(beta)) + l_ev_l(i)=lam(1) + l_ev_r(i)=lam(5) + xwind(i) = (lam(1) + lam(5)) / (abs(lam(1)) + abs(lam(5))) + xwind(i) = min(1.d0, max(-1.d0, xwind(i))) #define LEFTPLUS(x,xplus) xplus(i) = abs(xwind(i)) * LEFT1(x) + \ - (1.d0-abs(xwind(i))) * xplus(i) + (1.d0-abs(xwind(i))) * xplus(i) #define LEFTMINUS(x,xminus) xminus(i+1)= abs(xwind(i)) * LEFT1(x) + \ - (1.d0-abs(xwind(i))) * xminus(i+1) + (1.d0-abs(xwind(i))) * xminus(i+1) #define RIGHTPLUS(x,xplus) xplus(i) = abs(xwind(i)) * RIGHT1(x) + \ - (1.d0-abs(xwind(i))) * xplus(i) + (1.d0-abs(xwind(i))) * xplus(i) #define RIGHTMINUS(x,xminus) xminus(i+1)= abs(xwind(i)) * RIGHT1(x) + \ - (1.d0-abs(xwind(i))) * xminus(i+1) + (1.d0-abs(xwind(i))) * xminus(i+1) #define CHECK(x,xc) if (x(i+1) .gt. x(i)) then && xc=min(x(i+1),max(x(i),xc)) && else && xc=min(x(i),max(x(i+1),xc)) && endif -!!$ xwind(i)=0.d0 - if (xwind(i) .lt. 0.0d0) then - LEFTPLUS(rho, rhoplus) - LEFTMINUS(rho, rhominus) - LEFTPLUS(velx, velxplus) - LEFTMINUS(velx, velxminus) - LEFTPLUS(vely, velyplus) - LEFTMINUS(vely, velyminus) - LEFTPLUS(velz, velzplus) - LEFTMINUS(velz, velzminus) - if (poly .eq. 0) then - LEFTPLUS(eps, epsplus) - LEFTMINUS(eps, epsminus) - end if - else - RIGHTPLUS(rho, rhoplus) - RIGHTMINUS(rho, rhominus) - RIGHTPLUS(velx, velxplus) - RIGHTMINUS(velx, velxminus) - RIGHTPLUS(vely, velyplus) - RIGHTMINUS(vely, velyminus) - RIGHTPLUS(velz, velzplus) - RIGHTMINUS(velz, velzminus) - if (poly .eq. 0) then - RIGHTPLUS(eps, epsplus) - RIGHTMINUS(eps, epsminus) - end if - end if - CHECK(rho, rhoplus(i)) - CHECK(rho, rhominus(i+1)) - CHECK(velx, velxplus(i)) - CHECK(velx, velxminus(i+1)) - CHECK(vely, velyplus(i)) - CHECK(vely, velyminus(i+1)) - CHECK(velz, velzplus(i)) - CHECK(velz, velzminus(i+1)) - if (poly .eq. 0) then - CHECK(eps, epsplus(i)) - CHECK(eps, epsminus(i+1)) + !!$ xwind(i)=0.d0 + if (xwind(i) .lt. 0.0d0) then + LEFTPLUS(rho, rhoplus) + LEFTMINUS(rho, rhominus) + LEFTPLUS(velx, velxplus) + LEFTMINUS(velx, velxminus) + LEFTPLUS(vely, velyplus) + LEFTMINUS(vely, velyminus) + LEFTPLUS(velz, velzplus) + LEFTMINUS(velz, velzminus) + if (poly .eq. 0) then + LEFTPLUS(eps, epsplus) + LEFTMINUS(eps, epsminus) + end if + else + RIGHTPLUS(rho, rhoplus) + RIGHTMINUS(rho, rhominus) + RIGHTPLUS(velx, velxplus) + RIGHTMINUS(velx, velxminus) + RIGHTPLUS(vely, velyplus) + RIGHTMINUS(vely, velyminus) + RIGHTPLUS(velz, velzplus) + RIGHTMINUS(velz, velzminus) + if (poly .eq. 0) then + RIGHTPLUS(eps, epsplus) + RIGHTMINUS(eps, epsminus) + end if + end if + CHECK(rho, rhoplus(i)) + CHECK(rho, rhominus(i+1)) + CHECK(velx, velxplus(i)) + CHECK(velx, velxminus(i+1)) + CHECK(vely, velyplus(i)) + CHECK(vely, velyminus(i+1)) + CHECK(velz, velzplus(i)) + CHECK(velz, velzminus(i+1)) + if (poly .eq. 0) then + CHECK(eps, epsplus(i)) + CHECK(eps, epsminus(i+1)) + end if + !!$ if ((dir .eq. 1) .and. (ni .eq. 4) .and. (nj .eq. 4)) then + !!$ write (*,*) rhoplus(i), rhominus(i+1) + !!$ end if + end do + !!$ mppm debug output + if (ppm_mppm_debug_eigenvalues .gt. 0) then + if (dir .eq. 1) then + ev_l(:,ni,nj) = l_ev_l + ev_r(:,ni,nj) = l_ev_r + xw(:,ni,nj) = xwind + else if (dir .eq. 2) then + ev_l(ni,:,nj) = l_ev_l + ev_r(ni,:,nj) = l_ev_r + xw(ni,:,nj) = xwind + else if (dir .eq. 3) then + ev_l(ni,nj,:) = l_ev_l + ev_r(ni,nj,:) = l_ev_r + xw(ni,nj,:) = xwind + else + write (*,*) "flux direction not 1 to 3 ?" + end if end if -!!$ if ((dir .eq. 1) .and. (ni .eq. 4) .and. (nj .eq. 4)) then -!!$ write (*,*) rhoplus(i), rhominus(i+1) -!!$ end if - end do - !!$ mppm debug output - if (ppm_mppm_debug_eigenvalues .gt. 0) then - if (dir .eq. 1) then - ev_l(:,ni,nj) = l_ev_l - ev_r(:,ni,nj) = l_ev_r - xw(:,ni,nj) = xwind - else if (dir .eq. 2) then - ev_l(ni,:,nj) = l_ev_l - ev_r(ni,:,nj) = l_ev_r - xw(ni,:,nj) = xwind - else if (dir .eq. 3) then - ev_l(ni,nj,:) = l_ev_l - ev_r(ni,nj,:) = l_ev_r - xw(ni,nj,:) = xwind - else - write (*,*) "flux direction not 1 to 3 ?" - end if - end if + end if end if !!$ 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) @@ -320,54 +446,53 @@ trivial_rp = .true. end do - - 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) + 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) + 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 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 !!$ Monotonicity. See (1.10) of C&W. - -do i = GRHydro_stencil, nx - GRHydro_stencil + 1 #define MON(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&&\ @@ -393,16 +518,197 @@ do i = GRHydro_stencil, nx - GRHydro_stencil + 1 trivial_rp(i) = .false. &&\ end if - MON(rhominus,rho,rhoplus) - MON(velxminus,velx,velxplus) - MON(velyminus,vely,velyplus) - MON(velzminus,velz,velzplus) - end do - if (poly .eq. 0) then - do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - MON(epsminus,eps,epsplus) - end do - 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) \ + 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(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)) &&\ + 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 &&\ + else &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + 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 &&\ + 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! +#define MON_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(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)) &&\ + 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 &&\ + else &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + 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 &&\ + 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. + 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 + !! Constrain parabolic profiles, PPM 1984 + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + ! original Colella&Woodward monotonicity preservation + MON(rhominus,rho,rhoplus) + MON(velxminus,velx,velxplus) + MON(velyminus,vely,velyplus) + MON(velzminus,velz,velzplus) + end do + if (poly .eq. 0) then + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON(epsminus,eps,epsplus) + end do + end if + endif if (check_for_trivial_rp .eq. 0) then trivial_rp = .false. @@ -747,45 +1053,107 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & CCTK_REAL :: dpress2,w,flatten,dvel CCTK_REAL :: eta, etatilde -!!$ Average slopes delta_m(a). See (1.7) of Colella and Woodward, p.178 -!!$ This is the expression for an even grid. - + CCTK_REAL :: D2a, D2aL, D2aR, D2aC, D2aLim, rhi, threshold, rho_min, daplus, daminus, D3a, D3aLL, D3aL, D3aR, D2aLL, D2aRR, D3aMin, D3aMax + do i = 2, nx - 1 - dpress(i) = press(i+1) - press(i-1) + dpress(i) = press(i+1) - press(i-1) end do - do i = 2, nx - 1 - dY_e(i) = 0.5d0 * (Y_e(i+1) - Y_e(i-1)) -! d2Y_e(i,iY_e) = (Y_e(i+1) - 2.d0 * Y_e(i) + Y_e(i-1))! / 6.d0 / dx / dx -! ! since we use d2Y_e only for the condition d2Y_e(i+1)*d2Y_e(i-1)<0 -! ! the denominator is not necessary - enddo - - -!!$ Steepened slope. See (1.8) of Colella and Woodward, p.178 - - 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)) * & - 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))) - else - dmY_e(i) = 0.0d0 - endif - end do + + if (use_enhanced_PPM .eq. 0) then + + !!$ 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 + dY_e(i) = 0.5d0 * (Y_e(i+1) - Y_e(i-1)) + ! d2Y_e(i,iY_e) = (Y_e(i+1) - 2.d0 * Y_e(i) + Y_e(i-1))! / 6.d0 / dx / dx + ! ! since we use d2Y_e only for the condition d2Y_e(i+1)*d2Y_e(i-1)<0 + ! ! the denominator is not necessary + enddo + + + !!$ Steepened slope. See (1.8) of Colella and Woodward, p.178 + + 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)) * & + 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))) + else + dmY_e(i) = 0.0d0 + endif + end do + + !!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178 + + do i = 2, nx - 2 + Y_e_plus(i) = 0.5d0 * (Y_e(i) + Y_e(i+1)) + & + (dmY_e(i) - dmY_e(i+1)) / 6.d0 + Y_e_minus(i+1) = Y_e_plus(i) + enddo + else +!! This is the modified PPM algorithm by Colella & Sekora 2008. +!! 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, p.7071 +#undef APPROX_AT_CELL_INTERFACE_STENCIL4 +!!$ Initial boundary states (sixth order accurate). See (17) of Colella and Sekora 2008, p.7071 +#define APPROX_AT_CELL_INTERFACE_STENCIL4(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)) + +#undef APPROX_AT_CELL_INTERFACE + !!$ Initial boundary states (4th order accurate). See (16) of Colella and Sekora 2008, p.7071 +#define APPROX_AT_CELL_INTERFACE(a, ah) \ + ah = 7.0d0/12.0d0*(a(i)+a(i+1)) - 1.0d0/12.0d0*(a(i-1)+a(i+2)) + + +#undef LIMIT +#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 &&\ + 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)) &&\ + if (D2a*D2aR .ge. 0 .and. D2a*D2aL .ge. 0) 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)) &&\ + 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. + !! There, we also switch to lower order if we are near atmosphere values. + do i = 2, nx-2 + APPROX_AT_CELL_INTERFACE(Y_e, Y_e_plus(i)) + LIMIT(Y_e, Y_e_plus(i), enhanced_ppm_C2, Y_e_plus(i)) + Y_e_minus(i+1) = Y_e_plus(i) + end do + + else -!!$ Initial boundary states. See (1.9) of Colella and Woodward, p.178 + !! 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_STENCIL4(Y_e, Y_e_plus(i)) + LIMIT(Y_e, Y_e_plus(i), enhanced_ppm_C2, Y_e_plus(i)) + Y_e_minus(i+1) = Y_e_plus(i) + end do + endif + + end if - do i = 2, nx - 2 - Y_e_plus(i) = 0.5d0 * (Y_e(i) + Y_e(i+1)) + & - (dmY_e(i) - dmY_e(i+1)) / 6.d0 - Y_e_minus(i+1) = Y_e_plus(i) - enddo - !!$Discontinuity steepening. See (1.14-17) of C&W. !!$This is the detect routine which mat be activated with the ppm_detect parameter @@ -814,68 +1182,173 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & end if end do - if (PPM3) then - do i = 3, nx - 2 - flatten = tilde_flatten(i) - Y_e_plus(i) = flatten * Y_e_plus(i) & - + (1.d0 - flatten) * Y_e(i) - Y_e_minus(i) = flatten * Y_e_minus(i) & - + (1.d0 - flatten) * Y_e(i) - 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)) - Y_e_plus(i) = flatten * Y_e_plus(i) + & - (1.d0 - flatten) * Y_e(i) - Y_e_minus(i) = flatten * Y_e_minus(i) & - + (1.d0 - flatten) * Y_e(i) - end do + if (use_enhanced_ppm .eq. 0) then + if (PPM3) then + do i = 3, nx - 2 + flatten = tilde_flatten(i) + Y_e_plus(i) = flatten * Y_e_plus(i) & + + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + 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)) + Y_e_plus(i) = flatten * Y_e_plus(i) + & + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + end do + 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. +#define MON_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(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)) &&\ + 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 &&\ + else &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + 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 &&\ + 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! +#define MON_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(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)) &&\ + 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 &&\ + else &&\ + aplus(i) = a(i) + daplus * rhi &&\ + aminus(i) = a(i) - daminus * rhi &&\ + 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 &&\ + endif + + + if (use_enhanced_ppm .eq. 1) then + ! Constrain parabolic profiles + if (PPM3) then + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON_WITH_LOCAL_EXTREMUM(Y_e_minus,Y_e,Y_e_plus, enhanced_ppm_C2) + end do + else + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON_WITH_LOCAL_EXTREMUM_STENCIL4(Y_e_minus,Y_e,Y_e_plus, enhanced_ppm_C2, enhanced_ppm_C3) + end do + endif + + ! In PPM 2011/2008, flattening is applied after constraining parabolic profiles. + if (PPM3) then + do i = 3, nx - 2 + flatten = tilde_flatten(i) + Y_e_plus(i) = flatten * Y_e_plus(i) & + + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + 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)) + Y_e_plus(i) = flatten * Y_e_plus(i) + & + (1.d0 - flatten) * Y_e(i) + Y_e_minus(i) = flatten * Y_e_minus(i) & + + (1.d0 - flatten) * Y_e(i) + end do + end if + + else + !!$ Monotonicity. See (1.10) of C&W. + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + if (((Y_e_plus(i)-Y_e(i))* & + (Y_e(i)-Y_e_minus(i)) .le. 0.d0)) then + Y_e_minus(i) = Y_e(i) + Y_e_plus(i) = Y_e(i) + else if ((Y_e_plus(i) - Y_e_minus(i)) * (Y_e(i) - 0.5d0 * & + (Y_e_plus(i) + Y_e_minus(i))) > & + (Y_e_plus(i) - Y_e_minus(i))**2 / 6.d0) then + Y_e_minus(i) = 3.d0 * Y_e(i) - 2.0d0 * Y_e_plus(i) + else if ((Y_e_plus(i) - Y_e_minus(i)) * (Y_e(i) - 0.5d0 * & + (Y_e_plus(i) + Y_e_minus(i))) < & + -(Y_e_plus(i) - Y_e_minus(i))**2 / 6.0d0 ) then + Y_e_plus(i) = 3.d0 * Y_e(i) - 2.d0 * Y_e_minus(i) + end if + end do end if -!! Additional flattening a la Plewa & Mueller - do i = 2, nx - 1 - if ( ( Y_e(i+1) - Y_e(i) ) * & - ( Y_e(i) - Y_e(i-1) ) < 0.0d0 ) then - Y_eflat(i) = 1.0d0 - else - Y_eflat(i) = 0.0d0 - endif - enddo - - do i = 3, nx -2 - - Y_eflatomega = 0.5d0 * max(Y_eflat(i-1),2.0d0*Y_eflat(i), & - Y_eflat(i+1)) * ppm_omega_tracer - - Y_e_plus(i) = Y_eflatomega*Y_e(i) + & - (1.0d0 - Y_eflatomega)*Y_e_plus(i) - - Y_e_minus(i) = Y_eflatomega*Y_e(i) + & - (1.0d0 - Y_eflatomega)*Y_e_minus(i) - - enddo - - -!!$ Monotonicity. See (1.10) of C&W. - do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - if (((Y_e_plus(i)-Y_e(i))* & - (Y_e(i)-Y_e_minus(i)) .le. 0.d0)) then - Y_e_minus(i) = Y_e(i) - Y_e_plus(i) = Y_e(i) - else if ((Y_e_plus(i) - Y_e_minus(i)) * (Y_e(i) - 0.5d0 * & - (Y_e_plus(i) + Y_e_minus(i))) > & - (Y_e_plus(i) - Y_e_minus(i))**2 / 6.d0) then - Y_e_minus(i) = 3.d0 * Y_e(i) - 2.d0 * Y_e_plus(i) - else if ((Y_e_plus(i) - Y_e_minus(i)) * (Y_e(i) - 0.5d0 * & - (Y_e_plus(i) + Y_e_minus(i))) < & - -(Y_e_plus(i) - Y_e_minus(i))**2 / 6.d0 ) then - Y_e_plus(i) = 3.d0 * Y_e(i) - 2.d0 * Y_e_minus(i) - end if - end do - - - end subroutine SimplePPM_ye_1d - diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90 index 02ab359..b074bec 100644 --- a/src/GRHydro_PPMReconstruct_drv.F90 +++ b/src/GRHydro_PPMReconstruct_drv.F90 @@ -168,7 +168,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) w_lorentz(:,j,k), & 1, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_x_left, & GRHydro_mppm_eigenvalue_x_right, & - GRHydro_mppm_xwind) + GRHydro_mppm_xwind, GRHydro_rho_min) do i = 1, nx if (trivial_rp(i,j,k)) then SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bitsx, trivialx) @@ -198,7 +198,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) call SimplePPM_ye_1d(nx,CCTK_DELTA_SPACE(1),rho(:,j,k), & velx(:,j,k),vely(:,j,k),velz(:,j,k), & Y_e(:,j,k),Y_e_minus(:,j,k),Y_e_plus(:,j,k), & - press(:,j,k)) + press(:,j,k), GRHydro_rho_min) end do end do !$OMP END PARALLEL DO @@ -219,7 +219,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) w_lorentz(j,:,k), & 2, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_y_left, & GRHydro_mppm_eigenvalue_y_right, & - GRHydro_mppm_xwind) + GRHydro_mppm_xwind, GRHydro_rho_min) do i = 1, ny if (trivial_rp(j,i,k)) then SpaceMask_SetStateBitsF90(space_mask, j, i, k, type_bitsy, trivialy) @@ -249,7 +249,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) call SimplePPM_ye_1d(ny,CCTK_DELTA_SPACE(2),rho(j,:,k), & velx(j,:,k),vely(j,:,k),velz(j,:,k), & Y_e(j,:,k),Y_e_minus(j,:,k),Y_e_plus(j,:,k), & - press(j,:,k)) + press(j,:,k), GRHydro_rho_min) end do end do !$OMP END PARALLEL DO @@ -270,7 +270,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) w_lorentz(j,k,:), & 3, j, k, nx, ny, nz, GRHydro_mppm_eigenvalue_z_left, & GRHydro_mppm_eigenvalue_z_right, & - GRHydro_mppm_xwind) + GRHydro_mppm_xwind, GRHydro_rho_min) do i = 1, nz if (trivial_rp(j,k,i)) then SpaceMask_SetStateBitsF90(space_mask, j, k, i, type_bitsz, trivialz) @@ -301,7 +301,7 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) call SimplePPM_ye_1d(nz,CCTK_DELTA_SPACE(3),rho(j,k,:), & velx(j,k,:),vely(j,k,:),velz(j,k,:), & Y_e(j,k,:),Y_e_minus(j,k,:),Y_e_plus(j,k,:), & - press(j,k,:)) + press(j,k,:), GRHydro_rho_min) end do end do !$OMP END PARALLEL DO diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index cbcfd48..974d320 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -41,6 +41,20 @@ subroutine Reconstruction(CCTK_ARGUMENTS) CCTK_INT :: i,j,k CCTK_REAL :: local_min_tracer + ! save memory when MP is not used + ! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1 + TARGET lvel, vel + + ! save memory when MP is not used + CCTK_INT :: GRHydro_UseGeneralCoordinates + CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup + + if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then + vup => lvel + else + vup => vel + end if + if (CCTK_EQUALS(recon_method,"tvd")) then ! this handles MHD and non-MHD call GRHydro_TVDReconstruct_drv(CCTK_PASS_FTOF) @@ -79,12 +93,12 @@ subroutine Reconstruction(CCTK_ARGUMENTS) rhominus(i,j,k) = rho(i,j,k) epsplus(i,j,k) = eps(i,j,k) epsminus(i,j,k) = eps(i,j,k) - velxplus(i,j,k) = vel(i,j,k,1) - velyplus(i,j,k) = vel(i,j,k,2) - velzplus(i,j,k) = vel(i,j,k,3) - velxminus(i,j,k) = vel(i,j,k,1) - velyminus(i,j,k) = vel(i,j,k,2) - velzminus(i,j,k) = vel(i,j,k,3) + velxplus(i,j,k) = vup(i,j,k,1) + velyplus(i,j,k) = vup(i,j,k,2) + velzplus(i,j,k) = vup(i,j,k,3) + velxminus(i,j,k) = vup(i,j,k,1) + velyminus(i,j,k) = vup(i,j,k,2) + velzminus(i,j,k) = vup(i,j,k,3) end if if(evolve_tracer.ne.0) then where(tracerplus(i,j,k,:).le.local_min_tracer .or. & |