User: rhaas Date: 2013/01/14 08:23 AM
Modified: /trunk/src/ GRHydro_Con2PrimM.F90
Log: GRHydro: Call pointwise polytype upon con2primM failure. If it fails just use previous primitive values and get the cons from them. Basically freezing their values at previous time step.
From: Bruno Coutinho Mundim bcmsma@astro.rit.edu
File Changes:
Directory: /trunk/src/ ======================
File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +31 -19 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2013-01-14 14:23:19 UTC (rev 449) +++ trunk/src/GRHydro_Con2PrimM.F90 2013-01-14 14:23:21 UTC (rev 450) @@ -344,8 +344,8 @@ uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k))
- if(sdet.ge.sqrtdet_thr) then - if(GRHydro_C2P_failed(i,j,k).ne.0) then + if(GRHydro_C2P_failed(i,j,k).ne.0) then + if(sdet.ge.sqrtdet_thr) then GRHydro_C2P_failed(i,j,k) = 0
rho_tmp = rho(i,j,k) @@ -358,17 +358,38 @@ Bvecx_tmp = Bprim(i,j,k,1) Bvecy_tmp = Bprim(i,j,k,2) Bvecz_tmp = Bprim(i,j,k,3) + + xrho=1.0d0; xtemp=0.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_K = xpress(1);
- call GRHydro_Con2PrimM_pt(GRHydro_polytrope_handle, keytemp, & - GRHydro_eos_rf_prec, local_gam(1), dens(i,j,k), & - scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), & - Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3),xye(1), & - xtemp(1),rho_tmp,velx_tmp,vely_tmp,velz_tmp,& - eps_tmp,press_tmp,Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,& - w_lorentz_tmp,g11(i,j,k),g12(i,j,k),g13(i,j,k),& - g22(i,j,k),g23(i,j,k),g33(i,j,k), & + xrho=10.0d0; xeps=1.0d0 + call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& + xrho,xeps,xtemp,xye,xpress,keyerr,anyerr) + local_pgam=log(xpress(1)/local_K)/log(xrho(1)) + sc = local_K*dens(i,j,k) + + call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & + Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3),rho_tmp,& + velx_tmp,vely_tmp,velz_tmp,eps_tmp,press_tmp,& + Bvecx_tmp,Bvecy_tmp,Bvecz_tmp,b2,w_lorentz_tmp,& + g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k)) + + if(GRHydro_C2P_failed(i,j,k).ne.0) then + GRHydro_C2P_failed(i,j,k) = 0 + call prim2conM(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k), & + g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),det, & + dens(i,j,k),scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), & + tau(i,j,k),Bcons(i,j,k,1),Bcons(i,j,k,2),Bcons(i,j,k,3), & + rho(i,j,k),vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3), & + eps(i,j,k),press(i,j,k), & + Bprim(i,j,k,1),Bprim(i,j,k,2),Bprim(i,j,k,3),w_lorentz(i,j,k)) + cycle + end if end if end if
@@ -526,15 +547,6 @@ uxx,uxy,uxz,uyy,uyz,uzz,det, & epsnegative,GRHydro_C2P_failed(i,j,k))
- -! call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), & -! scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, & -! Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),& -! vup(i,j,k,1),vup(i,j,k,2),vup(i,j,k,3),eps(i,j,k),press(i,j,k),& -! Bprim(i,j,k,1), Bprim(i,j,k,2), Bprim(i,j,k,3),b2, w_lorentz(i,j,k),& -! g11(i,j,k),g12(i,j,k),g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k), & -! uxx,uxy,uxz,uyy,uyz,uzz,det, & - #endif
end if ! if (epsnegative .ne. 0) then