User: rhaas Date: 2012/07/17 12:08 PM
Modified: /trunk/src/ GRHydro_PPM.F90
Log: GRHydro: ePPM: Constrain reconstructed vel to below speed of light!
From: Christian Reisswig <reisswig@scriwalker.(none)>
File Changes:
Directory: /trunk/src/ ======================
File [modified]: GRHydro_PPM.F90 Delta lines: +23 -0 =================================================================== --- trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:22 UTC (rev 395) +++ trunk/src/GRHydro_PPM.F90 2012-07-17 17:08:25 UTC (rev 396) @@ -224,6 +224,7 @@ endif
+ if (PPM3) then
!! We initialize "plus" \equiv a_j+1/2 with (16) via APPROX_AT_CELL_INTERFACE, @@ -767,6 +768,24 @@
+!! Compute scalar product v^i v_i +#define VEL2(dvelx,dvely,dvelz) \ + (gxx(i)*dvelx*dvelx + gyy(i)*dvely*dvely + gzz(i) \ + *dvelz*dvelz + 2*gxy(i)*dvelx*dvely + 2*gxz(i)*dvelx *dvelz + 2*gyz(i) \ + *dvely*dvelz) + +!! Check if velocity is below speed of light. If not, reduce reconstruction to first order! +#define VELCHECK(vxminus,vx,vxplus,vyminus,vy,vyplus,vzminus,vz,vzplus) \ + if (1.0d0 .ge. VEL2(vxminus(i),vyminus(i),vzminus(i)) .or. 1.0d0 .ge. VEL2(vxplus(i),vyplus(i),vzplus(i))) then &&\ + vxminus(i) = vx(i) &&\ + vxplus(i) = vx(i) &&\ + vyminus(i) = vy(i) &&\ + vyplus(i) = vy(i) &&\ + vzminus(i) = vz(i) &&\ + vzplus(i) = vz(i) &&\ + endif + + if (use_enhanced_ppm .eq. 1) then !! Constrain parabolic profiles, PPM 2011/2008 if (PPM3) then @@ -775,6 +794,8 @@ 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) + ! Ensure that velocity is below speed of light! + VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus) end do if (poly .eq. 0) then do i = 3, nx - 2 @@ -788,6 +809,8 @@ 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) + ! Ensure that velocity is below speed of light! + VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus) end do if (poly .eq. 0) then do i = 4, nx - 3