diff options
author | knarf <knarf@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2010-04-06 19:15:52 +0000 |
---|---|---|
committer | knarf <knarf@c83d129a-5a75-4d5a-9c4d-ed3a5855bf45> | 2010-04-06 19:15:52 +0000 |
commit | 36b450e84061db36db617bae3b21636c09b1ac24 (patch) | |
tree | e734a4829c8ef1da412844fff0e44e5134656e0b | |
parent | d7d1a02dc4d3a3e1202185b4fffaf87bac2fcf32 (diff) |
unify code: make primitive2conservativegeneral() use prim2conpolytype()
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@96 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45
-rw-r--r-- | src/Whisky_Prim2Con.F90 | 95 |
1 files changed, 22 insertions, 73 deletions
diff --git a/src/Whisky_Prim2Con.F90 b/src/Whisky_Prim2Con.F90 index 2c823f5..4ab24f2 100644 --- a/src/Whisky_Prim2Con.F90 +++ b/src/Whisky_Prim2Con.F90 @@ -360,8 +360,10 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & w = 1.d0 / sqrt(1.d0 - w_tmp) endif - dpress = EOS_Pressure(handle, drho, deps) - deps = EOS_SpecificIntEnergy(handle, drho, dpress) + if (handle .ge. 0) then + dpress = EOS_Pressure(handle, drho, deps) + deps = EOS_SpecificIntEnergy(handle, drho, dpress) + end if vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz @@ -580,78 +582,25 @@ subroutine primitive2conservativegeneral(CCTK_ARGUMENTS) call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) - w_lorentzminus(i,j,k) = 1.d0 / sqrt(1.d0 - & - (gxxl * velxminus(i,j,k)**2 + & - gyyl * velyminus(i,j,k)**2 + & - gzzl * velzminus(i,j,k)**2 + & - 2.d0 * gxyl * velxminus(i,j,k) * velyminus(i,j,k) + & - 2.d0 * gxzl * velxminus(i,j,k) * velzminus(i,j,k) + & - 2.d0 * gyzl * velyminus(i,j,k) * velzminus(i,j,k) )) - vlowx = gxxl * velxminus(i,j,k) + & - gxyl * velyminus(i,j,k) + & - gxzl * velzminus(i,j,k) - vlowy = gxyl * velxminus(i,j,k) + & - gyyl * velyminus(i,j,k) + & - gyzl * velzminus(i,j,k) - vlowz = gxzl * velxminus(i,j,k) + & - gyzl * velyminus(i,j,k) + & - gzzl * velzminus(i,j,k) - densminus(i,j,k) = sqrt(avg_detl) * rhominus(i,j,k) * & - w_lorentzminus(i,j,k) - sxminus(i,j,k) = sqrt(avg_detl) * & - ( rhominus(i,j,k) * & - (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * & - w_lorentzminus(i,j,k)**2 * vlowx - syminus(i,j,k) = sqrt(avg_detl) * & - ( rhominus(i,j,k) * & - (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * & - w_lorentzminus(i,j,k)**2 * vlowy - szminus(i,j,k) = sqrt(avg_detl) * & - ( rhominus(i,j,k) * & - (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * & - w_lorentzminus(i,j,k)**2 * vlowz - tauminus(i,j,k) = sqrt(avg_detl) * & - ( ( rhominus(i,j,k) * & - (1.d0 + epsminus(i,j,k)) + pressminus(i,j,k) ) * & - w_lorentzminus(i,j,k)**2 - pressminus(i,j,k) ) - & - densminus(i,j,k) - - w_lorentzplus(i,j,k) = 1.d0 / sqrt(1.d0 - & - (gxxr * velxplus(i,j,k)**2 + & - gyyr * velyplus(i,j,k)**2 + & - gzzr * velzplus(i,j,k)**2 + & - 2.d0 * gxyr * velxplus(i,j,k) * velyplus(i,j,k) + & - 2.d0 * gxzr * velxplus(i,j,k) * velzplus(i,j,k) + & - 2.d0 * gyzr * velyplus(i,j,k) * velzplus(i,j,k) )) - vlowx = gxxr * velxplus(i,j,k) + & - gxyr * velyplus(i,j,k) + & - gxzr * velzplus(i,j,k) - vlowy = gxyr * velxplus(i,j,k) + & - gyyr * velyplus(i,j,k) + & - gyzr * velzplus(i,j,k) - vlowz = gxzr * velxplus(i,j,k) + & - gyzr * velyplus(i,j,k) + & - gzzr * velzplus(i,j,k) - densplus(i,j,k) = sqrt(avg_detr) * rhoplus(i,j,k) * & - w_lorentzplus(i,j,k) - sxplus(i,j,k) = sqrt(avg_detr) * & - ( rhoplus(i,j,k) * & - (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * & - w_lorentzplus(i,j,k)**2 * vlowx - syplus(i,j,k) = sqrt(avg_detr) * & - ( rhoplus(i,j,k) * & - (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * & - w_lorentzplus(i,j,k)**2 * vlowy - szplus(i,j,k) = sqrt(avg_detr) * & - ( rhoplus(i,j,k) * & - (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * & - w_lorentzplus(i,j,k)**2 * vlowz - tauplus(i,j,k) = sqrt(avg_detr) * & - ( ( rhoplus(i,j,k) * & - (1.d0 + epsplus(i,j,k)) + pressplus(i,j,k) ) * & - w_lorentzplus(i,j,k)**2 - pressplus(i,j,k) ) - & - densplus(i,j,k) + call prim2conpolytype(-1, gxxl,gxyl,gxzl,& + gyyl,gyzl,gzzl, & + avg_detl,densminus(i,j,k),sxminus(i,j,k),& + syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),rhominus(i,j,k), & + velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& + epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i, j, k)) + call prim2conpolytype(-1, gxxr,gxyr,gxzr,& + gyyr,gyzr,gzzr, & + avg_detr, densplus(i,j,k),sxplus(i,j,k),& + syplus(i,j,k),szplus(i,j ,k),tauplus(i,j,k),& + rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& + velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& + w_lorentzplus(i,j,k)) + if (densminus(i,j,k) .ne. densminus(i,j,k)) then + !$OMP CRITICAL + call CCTK_WARN(0, "NaN in densminus(i,j,k) (Prim2Con)") + !$OMP END CRITICAL + endif end do end do end do |