diff options
Diffstat (limited to 'src/GRHydro_PPM.F90')
| -rw-r--r-- | src/GRHydro_PPM.F90 | 1646 |
1 files changed, 823 insertions, 823 deletions
diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index a6c46dc..f71a0e1 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -130,13 +130,13 @@ trivial_rp = .true. #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 ( (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 @@ -146,57 +146,57 @@ trivial_rp = .true. !!$ 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 * (vx(i+1) - vx(i-1)) - dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1)) - dvelz(i) = 0.5d0 * (vz(i+1) - vz(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 + drho(i) = 0.5d0 * (rho(i+1) - rho(i-1)) + dvelx(i) = 0.5d0 * (vx(i+1) - vx(i-1)) + dvely(i) = 0.5d0 * (vy(i+1) - vy(i-1)) + dvelz(i) = 0.5d0 * (vz(i+1) - vz(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 + 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 do i = 2, nx - 1 - STEEP(rho, drho, dmrho) - STEEP(vx, dvelx, dmvelx) - STEEP(vy, dvely, dmvely) - STEEP(vz, dvelz, dmvelz) + STEEP(rho, drho, dmrho) + STEEP(vx, dvelx, dmvelx) + STEEP(vy, dvely, dmvely) + STEEP(vz, dvelz, dmvelz) end do if (poly .eq. 0) then - do i = 2, nx - 1 - STEEP(eps, deps, dmeps) - end do + 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 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 * (vx(i) + vx(i+1)) + & - (dmvelx(i) - dmvelx(i+1)) / 6.d0 - velxminus(i+1) = velxplus(i) - velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + & - (dmvely(i) - dmvely(i+1)) / 6.d0 - velyminus(i+1) = velyplus(i) - velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + & - (dmvelz(i) - dmvelz(i+1)) / 6.d0 - velzminus(i+1) = velzplus(i) + 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 * (vx(i) + vx(i+1)) + & + (dmvelx(i) - dmvelx(i+1)) / 6.d0 + velxminus(i+1) = velxplus(i) + velyplus(i) = 0.5d0 * (vy(i) + vy(i+1)) + & + (dmvely(i) - dmvely(i+1)) / 6.d0 + velyminus(i+1) = velyplus(i) + velzplus(i) = 0.5d0 * (vz(i) + vz(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 + 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. @@ -213,32 +213,32 @@ trivial_rp = .true. #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 &&\ + 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 &&\ + 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 LIMIT_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 &&\ + 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 &&\ + 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 @@ -266,11 +266,11 @@ trivial_rp = .true. 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(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) - epsminus(i+1) = epsplus(i) - end do + do i = 2, nx-2 + APPROX_AT_CELL_INTERFACE(eps, epsplus(i)) + LIMIT_EPS(eps, epsplus(i), enhanced_ppm_C2, epsplus(i)) + epsminus(i+1) = epsplus(i) + end do end if else @@ -295,11 +295,11 @@ trivial_rp = .true. 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 + 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 @@ -323,36 +323,36 @@ trivial_rp = .true. 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 + ! 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 + 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 @@ -366,95 +366,95 @@ trivial_rp = .true. 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))) + 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)) - end if + 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 + 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 @@ -481,46 +481,46 @@ trivial_rp = .true. 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) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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 + 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) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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(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. - 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) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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 + do i = 4, nx - 3 + 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. + 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) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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 @@ -557,50 +557,50 @@ trivial_rp = .true. 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 &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 &&\ + 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*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ endif !! Monotonicity of PPM 2011 (McCorquodale & Colella 2011), Sec. 2.4.1, Eq. 23-34. @@ -611,69 +611,69 @@ trivial_rp = .true. 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 &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ 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 &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 @@ -684,40 +684,40 @@ trivial_rp = .true. 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 &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 &&\ + 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*daminus &&\ - else if (abs(daminus) .ge. 2.0d0*abs(daplus)) then &&\ - aminus(i) = a(i) - 2.0d0*daplus &&\ - end if &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ endif @@ -729,133 +729,133 @@ trivial_rp = .true. 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 &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 &&\ + trivial_rp(i-1) = .false. &&\ + trivial_rp(i) = .false. &&\ 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 &&\ - trivial_rp(i-1) = .false. &&\ - trivial_rp(i) = .false. &&\ + 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 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,vx,velxplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velyminus,vy,velyplus, enhanced_ppm_C2) - MON_WITH_LOCAL_EXTREMUM(velzminus,vz,velzplus, enhanced_ppm_C2) - end do - if (poly .eq. 0) then - do i = 3, nx - 2 - MON_WITH_LOCAL_EXTREMUM_EPS(epsminus,eps,epsplus, enhanced_ppm_C2) - end do - end if - + do i = 3, nx - 2 + MON_WITH_LOCAL_EXTREMUM(rhominus,rho,rhoplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velxminus,vx,velxplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velyminus,vy,velyplus, enhanced_ppm_C2) + MON_WITH_LOCAL_EXTREMUM(velzminus,vz,velzplus, enhanced_ppm_C2) + end do + if (poly .eq. 0) then + do i = 3, nx - 2 + MON_WITH_LOCAL_EXTREMUM_EPS(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,vx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vy,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) - MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,vz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) - end do - if (poly .eq. 0) then - do i = 4, nx - 3 - MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(epsminus,eps,epsplus, enhanced_ppm_C2, enhanced_ppm_C3) - end do - endif + 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,vx,velxplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velyminus,vy,velyplus, enhanced_ppm_C2, enhanced_ppm_C3) + MON_WITH_LOCAL_EXTREMUM_STENCIL4(velzminus,vz,velzplus, enhanced_ppm_C2, enhanced_ppm_C3) + end do + if (poly .eq. 0) then + do i = 4, nx - 3 + MON_WITH_LOCAL_EXTREMUM_STENCIL4_EPS(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 - 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) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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 + 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) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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(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. - 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) * vx(i) - velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) - velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) - velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) - velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) - velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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 + do i = 4, nx - 3 + 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. + 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) * vx(i) + velxminus(i) = flatten * velxminus(i) + (1.d0 - flatten) * vx(i) + velyplus(i) = flatten * velyplus(i) + (1.d0 - flatten) * vy(i) + velyminus(i) = flatten * velyminus(i) + (1.d0 - flatten) * vy(i) + velzplus(i) = flatten * velzplus(i) + (1.d0 - flatten) * vz(i) + velzminus(i) = flatten * velzminus(i) + (1.d0 - flatten) * vz(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 @@ -863,16 +863,16 @@ trivial_rp = .true. 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,vx,velxplus) - MON(velyminus,vy,velyplus) - MON(velzminus,vz,velzplus) + ! original Colella&Woodward monotonicity preservation + MON(rhominus,rho,rhoplus) + MON(velxminus,vx,velxplus) + MON(velyminus,vy,velyplus) + MON(velzminus,vz,velzplus) end do if (poly .eq. 0) then - do i = GRHydro_stencil, nx - GRHydro_stencil + 1 - MON(epsminus,eps,epsplus) - end do + do i = GRHydro_stencil, nx - GRHydro_stencil + 1 + MON(epsminus,eps,epsplus) + end do end if endif @@ -1067,13 +1067,13 @@ subroutine SimplePPM_temperature_1d(& #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 ( (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 @@ -1117,32 +1117,32 @@ subroutine SimplePPM_temperature_1d(& #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 &&\ + 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 &&\ + 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 &&\ + 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&&\ + 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 @@ -1202,18 +1202,18 @@ subroutine SimplePPM_temperature_1d(& 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) + 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 + 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)) + 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 do end if endif @@ -1240,46 +1240,46 @@ subroutine SimplePPM_temperature_1d(& 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 &&\ + 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 &&\ + 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. @@ -1290,65 +1290,65 @@ subroutine SimplePPM_temperature_1d(& 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 &&\ + 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 &&\ + 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 @@ -1358,36 +1358,36 @@ subroutine SimplePPM_temperature_1d(& 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 &&\ + 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 &&\ + 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 @@ -1399,84 +1399,84 @@ subroutine SimplePPM_temperature_1d(& 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 &&\ + 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 &&\ + 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 + do i = 3, nx - 2 MONT_WITH_LOCAL_EXTREMUM_EPS(tempminus,temperature,tempplus, enhanced_ppm_C2) enddo else - do i = 4, nx - 3 + 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) + 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 + 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)) + 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 do end if endif @@ -1484,7 +1484,7 @@ subroutine SimplePPM_temperature_1d(& 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 + ! original Colella&Woodward monotonicity preservation MONT(tempminus,temperature,tempplus) enddo endif @@ -1850,38 +1850,38 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & #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 &&\ + 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(one, D2a)*min(C*abs(D2aL), C*abs(D2aR), abs(D2a))/3.0d0 &&\ - else &&\ - alim = 0.5d0*(a(i)+a(i+1)) &&\ - end if &&\ + 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(one, 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 - + !! 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 - !! 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 + !! 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 @@ -1917,22 +1917,22 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & 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 + 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(one, -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 + do i = 4, nx - 3 + s=sign(one, -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 @@ -1944,46 +1944,46 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & 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 &&\ + 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 &&\ + 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 @@ -1997,36 +1997,36 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & 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 &&\ + 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*(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) .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 @@ -2036,33 +2036,33 @@ subroutine SimplePPM_ye_1d(nx,dx,rho,velx,vely,velz, & 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 + 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 + 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 + 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(one, -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 + do i = 4, nx - 3 + s=sign(one, -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 |