[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 416)
rhaas at tapir.caltech.edu
rhaas at tapir.caltech.edu
Mon Aug 27 14:19:38 CDT 2012
User: rhaas
Date: 2012/08/27 02:19 PM
Modified:
/trunk/
param.ccl
/trunk/src/
GRHydro_PPM.F90, GRHydro_Reconstruct.F90
Log:
GRHydro: Option to reconstruct W*vel (W = Lorentz factor) instead of just the primitive velocity 'vel'.
This ensures that the Lorentz factor can never become unphysical during reconstruction.
This patch is necessary for NS collapse when ePPM+Refluxing is used.
The idea is due to Roland.
From: Christian Reisswig <reisswig at scriwalker.(none)>
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_PPM.F90
Delta lines: +147 -123
===================================================================
--- trunk/src/GRHydro_PPM.F90 2012-08-27 19:19:35 UTC (rev 415)
+++ trunk/src/GRHydro_PPM.F90 2012-08-27 19:19:38 UTC (rev 416)
@@ -106,7 +106,25 @@
logical :: cond
+ CCTK_REAL, dimension(nx) :: vx,vy,vz
+
+
+ !!$ reconstruct w^i = w_lorentz*v^i to ensure slower than light
+ !!$ speeds?
+ if (reconstruct_Wv.ne.0) then
+ ! all velocity like quantities are now w_lorentz*velocity. We will convert
+ ! back to ordinary velocity at the end, using w_lorentz = sqrt(1 + g_{ij}
+ ! w^i w^j).
+ vx = w_lorentz*velx
+ vy = w_lorentz*vely
+ vz = w_lorentz*velz
+ else
+ vx = velx
+ vy = vely
+ vz = velz
+ end if
+
!!$ Initially, all the Riemann problems will be trivial
trivial_rp = .true.
@@ -129,9 +147,9 @@
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))
+ 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
@@ -147,9 +165,9 @@
do i = 2, nx - 1
STEEP(rho, drho, dmrho)
- STEEP(velx, dvelx, dmvelx)
- STEEP(vely, dvely, dmvely)
- STEEP(velz, dvelz, dmvelz)
+ STEEP(vx, dvelx, dmvelx)
+ STEEP(vy, dvely, dmvely)
+ STEEP(vz, dvelz, dmvelz)
end do
if (poly .eq. 0) then
do i = 2, nx - 1
@@ -163,13 +181,13 @@
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)) + &
+ 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 * (vely(i) + vely(i+1)) + &
+ 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 * (velz(i) + velz(i+1)) + &
+ velzplus(i) = 0.5d0 * (vz(i) + vz(i+1)) + &
(dmvelz(i) - dmvelz(i+1)) / 6.d0
velzminus(i+1) = velzplus(i)
end do
@@ -235,16 +253,16 @@
LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i))
rhominus(i+1) = rhoplus(i)
- APPROX_AT_CELL_INTERFACE(velx, velxplus(i))
- LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i))
+ APPROX_AT_CELL_INTERFACE(vx, velxplus(i))
+ LIMIT(vx, velxplus(i), enhanced_ppm_C2, velxplus(i))
velxminus(i+1) = velxplus(i)
- APPROX_AT_CELL_INTERFACE(vely, velyplus(i))
- LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i))
+ APPROX_AT_CELL_INTERFACE(vy, velyplus(i))
+ LIMIT(vy, velyplus(i), enhanced_ppm_C2, velyplus(i))
velyminus(i+1) = velyplus(i)
- APPROX_AT_CELL_INTERFACE(velz, velzplus(i))
- LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i))
+ APPROX_AT_CELL_INTERFACE(vz, velzplus(i))
+ LIMIT(vz, velzplus(i), enhanced_ppm_C2, velzplus(i))
velzminus(i+1) = velzplus(i)
end do
if (poly .eq. 0) then
@@ -264,16 +282,16 @@
LIMIT(rho, rhoplus(i), enhanced_ppm_C2, rhoplus(i))
rhominus(i+1) = rhoplus(i)
- APPROX_AT_CELL_INTERFACE_STENCIL4(velx, velxplus(i))
- LIMIT(velx, velxplus(i), enhanced_ppm_C2, velxplus(i))
+ APPROX_AT_CELL_INTERFACE_STENCIL4(vx, velxplus(i))
+ LIMIT(vx, velxplus(i), enhanced_ppm_C2, velxplus(i))
velxminus(i+1) = velxplus(i)
- APPROX_AT_CELL_INTERFACE_STENCIL4(vely, velyplus(i))
- LIMIT(vely, velyplus(i), enhanced_ppm_C2, velyplus(i))
+ APPROX_AT_CELL_INTERFACE_STENCIL4(vy, velyplus(i))
+ LIMIT(vy, velyplus(i), enhanced_ppm_C2, velyplus(i))
velyminus(i+1) = velyplus(i)
- APPROX_AT_CELL_INTERFACE_STENCIL4(velz, velzplus(i))
- LIMIT(velz, velzplus(i), enhanced_ppm_C2, velzplus(i))
+ APPROX_AT_CELL_INTERFACE_STENCIL4(vz, velzplus(i))
+ LIMIT(vz, velzplus(i), enhanced_ppm_C2, velzplus(i))
velzminus(i+1) = velzplus(i)
end do
if (poly .eq. 0) then
@@ -444,7 +462,7 @@
!!$ 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)
+ dvel = vx(i+1) - vx(i-1)
if ( (abs(dpress(i)) > ppm_epsilon * min(press(i-1),press(i+1))) .and. &
(dvel < 0.d0) ) then
w = 1.d0
@@ -471,12 +489,12 @@
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)
+ 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)
@@ -492,12 +510,12 @@
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)
+ 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)
@@ -767,50 +785,14 @@
endif
-! Reconstruct quantity on plus face
-#define AVGP(a) \
- 0.5*(a(i)+a(i+1))
-
-! Reconstruct quantity on minus face
-#define AVGM(a) \
- 0.5*(a(i)+a(i-1))
-
-
-!! Compute scalar product v^i v_i on plus face
-#define VEL2P(dvelx,dvely,dvelz) \
- (AVGP(gxx)*dvelx*dvelx + AVGP(gyy)*dvely*dvely + AVGP(gzz) \
- *dvelz*dvelz + 2*AVGP(gxy)*dvelx*dvely + 2*AVGP(gxz)*dvelx *dvelz + 2*AVGP(gyz) \
- *dvely*dvelz)
-
-!! Compute scalar product v^i v_i on minus face
-#define VEL2M(dvelx,dvely,dvelz) \
- (AVGM(gxx)*dvelx*dvelx + AVGM(gyy)*dvely*dvely + AVGM(gzz) \
- *dvelz*dvelz + 2*AVGM(gxy)*dvelx*dvely + 2*AVGM(gxz)*dvelx *dvelz + 2*AVGM(gyz) \
- *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 (VEL2M(vxminus(i),vyminus(i),vzminus(i)) .ge. (1.0d0-enhanced_ppm_C4) .or. VEL2P(vxplus(i),vyplus(i),vzplus(i)) .ge. (1.0d0-enhanced_ppm_C4)) 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
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)
- ! Ensure that velocity is below speed of light!
- VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus)
+ 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
@@ -821,11 +803,9 @@
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)
- ! Ensure that velocity is below speed of light!
- VELCHECK(velxminus,velx,velxplus,velyminus,vely,velyplus,velzminus,velz,velzplus)
+ 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
@@ -844,12 +824,12 @@
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)
+ 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)
@@ -865,12 +845,12 @@
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)
+ 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)
@@ -885,9 +865,9 @@
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)
+ 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
@@ -919,20 +899,20 @@
if (cond) then
rhominus(i)=rho(i)
rhoplus(i)=rho(i)
- velxminus(i)=velx(i)
- velxplus(i)=velx(i)
- velyminus(i)=vely(i)
- velyplus(i)=vely(i)
- velzminus(i)=velz(i)
- velzplus(i)=velz(i)
+ velxminus(i)=vx(i)
+ velxplus(i)=vx(i)
+ velyminus(i)=vy(i)
+ velyplus(i)=vy(i)
+ velzminus(i)=vz(i)
+ velzplus(i)=vz(i)
rhominus(i-1)=rho(i)
rhoplus(i-1)=rho(i)
- velxminus(i-1)=velx(i)
- velxplus(i-1)=velx(i)
- velyminus(i-1)=vely(i)
- velyplus(i-1)=vely(i)
- velzminus(i-1)=velz(i)
- velzplus(i-1)=velz(i)
+ velxminus(i-1)=vx(i)
+ velxplus(i-1)=vx(i)
+ velyminus(i-1)=vy(i)
+ velyplus(i-1)=vy(i)
+ velzminus(i-1)=vz(i)
+ velzplus(i-1)=vz(i)
if (poly .eq. 0) then
epsminus(i)=eps(i)
epsplus(i)=eps(i)
@@ -946,9 +926,9 @@
end if
if (cond) then
call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
- call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i))
- call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i))
- call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i))
+ call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i))
+ call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i))
+ call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i))
if (poly .eq. 0) then
call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
end if
@@ -961,20 +941,20 @@
if (cond) then
rhominus(i)=rho(i)
rhoplus(i)=rho(i)
- velxminus(i)=velx(i)
- velxplus(i)=velx(i)
- velyminus(i)=vely(i)
- velyplus(i)=vely(i)
- velzminus(i)=velz(i)
- velzplus(i)=velz(i)
+ velxminus(i)=vx(i)
+ velxplus(i)=vx(i)
+ velyminus(i)=vy(i)
+ velyplus(i)=vy(i)
+ velzminus(i)=vz(i)
+ velzplus(i)=vz(i)
rhominus(i+1)=rho(i)
rhoplus(i+1)=rho(i)
- velxminus(i+1)=velx(i)
- velxplus(i+1)=velx(i)
- velyminus(i+1)=vely(i)
- velyplus(i+1)=vely(i)
- velzminus(i+1)=velz(i)
- velzplus(i+1)=velz(i)
+ velxminus(i+1)=vx(i)
+ velxplus(i+1)=vx(i)
+ velyminus(i+1)=vy(i)
+ velyplus(i+1)=vy(i)
+ velzminus(i+1)=vz(i)
+ velzplus(i+1)=vz(i)
if (poly .eq. 0) then
epsminus(i)=eps(i)
epsplus(i)=eps(i)
@@ -988,9 +968,9 @@
end if
if (cond) then
call PPM_TVD(rho(i-1), rho(i), rho(i+1), rhominus(i), rhoplus(i))
- call PPM_TVD(velx(i-1), velx(i), velx(i+1), velxminus(i), velxplus(i))
- call PPM_TVD(vely(i-1), vely(i), vely(i+1), velyminus(i), velyplus(i))
- call PPM_TVD(velz(i-1), velz(i), velz(i+1), velzminus(i), velzplus(i))
+ call PPM_TVD(vx(i-1), vx(i), vx(i+1), velxminus(i), velxplus(i))
+ call PPM_TVD(vy(i-1), vy(i), vy(i+1), velyminus(i), velyplus(i))
+ call PPM_TVD(vz(i-1), vz(i), vz(i+1), velzminus(i), velzplus(i))
if (poly .eq. 0) then
call PPM_TVD(eps(i-1), eps(i), eps(i+1), epsminus(i), epsplus(i))
end if
@@ -998,6 +978,50 @@
end if
end if
end do
+
+ !!$ transform back to vel if Wv was used!
+ if (reconstruct_Wv.ne.0) then
+ do i = grhydro_stencil, nx - grhydro_stencil + 1
+ ! divide out the Loretnz factor obtained from w_lorentz =
+ ! sqrt(1+g_{ij} w^i w^j) for both the
+ ! plus and minus quantities this should by construction ensure
+ ! that any Lorentz factor calculated from them later on is
+ ! physical (ie. > 1.d0)
+ 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) )
+ w = sqrt( 1.d0 + agxx*velxminus(i)*velxminus(i) + &
+ agyy*velyminus(i)*velyminus(i) &
+ + agzz*velzminus(i)*velzminus(i) + &
+ 2.d0*agxy*velxminus(i)*velyminus(i) &
+ + 2.d0*agxz*velxminus(i)*velzminus(i) + &
+ 2.d0*agyz*velyminus(i)*velzminus(i) )
+ velxminus(i) = velxminus(i)/w
+ velyminus(i) = velyminus(i)/w
+ velzminus(i) = velzminus(i)/w
+
+ 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) )
+ w = sqrt( 1.d0 + agxx*velxplus(i)*velxplus(i) + &
+ agyy*velyplus(i)*velyplus(i) &
+ + agzz*velzplus(i)*velzplus(i) + &
+ 2.d0*agxy*velxplus(i)*velyplus(i) &
+ + 2.d0*agxz*velxplus(i)*velzplus(i) + &
+ 2.d0*agyz*velyplus(i)*velzplus(i) )
+ velxplus(i) = velxplus(i)/w
+ velyplus(i) = velyplus(i)/w
+ velzplus(i) = velzplus(i)/w
+ end do
+ end if
+
+
return
end subroutine SimplePPM_1d
File [modified]: GRHydro_Reconstruct.F90
Delta lines: +68 -7
===================================================================
--- trunk/src/GRHydro_Reconstruct.F90 2012-08-27 19:19:35 UTC (rev 415)
+++ trunk/src/GRHydro_Reconstruct.F90 2012-08-27 19:19:38 UTC (rev 416)
@@ -36,6 +36,8 @@
! save memory when MP is not used
! TARGET as to be before DECLARE_CCTK_ARGUMENTS for gcc 4.1
+ TARGET gaa, gab, gac, gbb, gbc, gcc
+ TARGET gxx, gxy, gxz, gyy, gyz, gzz
TARGET lvel, vel
DECLARE_CCTK_ARGUMENTS
@@ -45,14 +47,28 @@
CCTK_INT :: i,j,k
CCTK_REAL :: local_min_tracer
CCTK_INT :: type_bits, not_trivial
+ CCTK_REAL :: agxx, agxy, agxz, agyy, agyz, agzz, w
! save memory when MP is not used
CCTK_INT :: GRHydro_UseGeneralCoordinates
+ CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33
CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup
if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
+ g11 => gaa
+ g12 => gab
+ g13 => gac
+ g22 => gbb
+ g23 => gbc
+ g33 => gcc
vup => lvel
else
+ g11 => gxx
+ g12 => gxy
+ g13 => gxz
+ g22 => gyy
+ g23 => gyz
+ g33 => gzz
vup => vel
end if
@@ -101,7 +117,7 @@
call CCTK_WARN(0, "Flux direction not x,y,z")
end if
- !$OMP PARALLEL DO PRIVATE(i,j,k)
+ !$OMP PARALLEL DO PRIVATE(i,j,k, agxx, agxy, agxz, agyy, agyz, agzz, w)
do k = GRHydro_stencil, cctk_lsh(3) - GRHydro_stencil + 1
do j = GRHydro_stencil, cctk_lsh(2) - GRHydro_stencil + 1
do i = GRHydro_stencil, cctk_lsh(1) - GRHydro_stencil + 1
@@ -112,12 +128,57 @@
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) = 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)
+ if (reconstruct_Wv.ne.0) then
+ ! divide out the Loretnz factor for both the
+ ! plus and minus quantities this should by construction ensure
+ ! that any Lorentz factor calculated from them later on is
+ ! physical (ie. > 1.d0)
+ velxplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1)
+ velyplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2)
+ velzplus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3)
+ agxx = 0.5d0 * (g11(i,j,k) + g11(i+xoffset,j+yoffset,k+zoffset))
+ agxy = 0.5d0 * (g12(i,j,k) + g12(i+xoffset,j+yoffset,k+zoffset))
+ agxz = 0.5d0 * (g13(i,j,k) + g13(i+xoffset,j+yoffset,k+zoffset))
+ agyy = 0.5d0 * (g22(i,j,k) + g22(i+xoffset,j+yoffset,k+zoffset))
+ agyz = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset))
+ agzz = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset))
+ w = sqrt( 1.d0 + agxx*velxplus(i,j,k)*velxplus(i,j,k) + &
+ agyy*velyplus(i,j,k)*velyplus(i,j,k) &
+ + agzz*velzplus(i,j,k)*velzplus(i,j,k) + &
+ 2.d0*agxy*velxplus(i,j,k)*velyplus(i,j,k) &
+ + 2.d0*agxz*velxplus(i,j,k)*velzplus(i,j,k) + &
+ 2.d0*agyz*velyplus(i,j,k)*velzplus(i,j,k) )
+ velxplus(i,j,k) = velxplus(i,j,k)/w
+ velyplus(i,j,k) = velyplus(i,j,k)/w
+ velzplus(i,j,k) = velzplus(i,j,k)/w
+
+ velxminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,1)
+ velyminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,2)
+ velzminus(i,j,k) = w_lorentz(i,j,k)*vup(i,j,k,3)
+ agxx = 0.5d0 * (g11(i,j,k) + g11(i-xoffset,j-yoffset,k-zoffset))
+ agxy = 0.5d0 * (g12(i,j,k) + g12(i-xoffset,j-yoffset,k-zoffset))
+ agxz = 0.5d0 * (g13(i,j,k) + g13(i-xoffset,j-yoffset,k-zoffset))
+ agyy = 0.5d0 * (g22(i,j,k) + g22(i-xoffset,j-yoffset,k-zoffset))
+ agyz = 0.5d0 * (g23(i,j,k) + g23(i-xoffset,j-yoffset,k-zoffset))
+ agzz = 0.5d0 * (g33(i,j,k) + g33(i-xoffset,j-yoffset,k-zoffset))
+ w = sqrt( 1.d0 + agxx*velxminus(i,j,k)*velxminus(i,j,k) + &
+ agyy*velyminus(i,j,k)*velyminus(i,j,k) &
+ + agzz*velzminus(i,j,k)*velzminus(i,j,k) + &
+ 2.d0*agxy*velxminus(i,j,k)*velyminus(i,j,k) &
+ + 2.d0*agxz*velxminus(i,j,k)*velzminus(i,j,k) + &
+ 2.d0*agyz*velyminus(i,j,k)*velzminus(i,j,k) )
+ velxminus(i,j,k) = velxminus(i,j,k)/w
+ velyminus(i,j,k) = velyminus(i,j,k)/w
+ velzminus(i,j,k) = velzminus(i,j,k)/w
+ else
+ ! This is the standard way of doing it
+ 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)
+ endif
if(evolve_y_e.ne.0) then
y_e_plus(i,j,k) = y_e(i,j,k)
y_e_minus(i,j,k) = y_e(i,j,k)
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +3 -0
===================================================================
--- trunk/param.ccl 2012-08-27 19:19:35 UTC (rev 415)
+++ trunk/param.ccl 2012-08-27 19:19:38 UTC (rev 416)
@@ -231,6 +231,9 @@
0:1 :: "must be greater than or equal to 0, and not larger than 1"
} 0.0
+boolean reconstruct_Wv "Reconstruct the primitive velocity W_Lorentz*vel, rather than just vel." STEERABLE=ALWAYS
+{
+} no
int eno_order "The order of accuracy of the ENO reconstruction"
{
More information about the Commits
mailing list