From 76ba623bb876010bae337587f20a440531956fe0 Mon Sep 17 00:00:00 2001 From: rhaas Date: Thu, 15 Sep 2011 16:43:29 +0000 Subject: add more workarounds in Con2Prim for hot EOS * fix con2prim issues (pertaining to error messages) * fix OMP CRITICAL issues in Prim2Con and Con2Prim original commits by Christian Ott (cott) git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@267 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Con2Prim.F90 | 133 +++++++++++++++++++++++++++++++---------------- src/GRHydro_Prim2Con.F90 | 12 ++--- 2 files changed, 94 insertions(+), 51 deletions(-) (limited to 'src') diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 2ba43cd..8279057 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -52,7 +52,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) integer :: i, j, k, itracer, nx, ny, nz CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, pmin, epsmin logical :: epsnegative - character(len=100) warnline + character*256 :: warnline CCTK_INT :: type_bits, atmosphere CCTK_INT :: type2_bits @@ -88,7 +88,8 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) GRHydro_rho_min,xeps,xtemp,xye,pmin,epsmin,keyerr,anyerr) !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,& - !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp) + !$OMP uxx, uxy, uxz, uyy, uyz, uzz, det, epsnegative, anyerr, keyerr, keytemp,& + !$OMP warnline) do k = 1, nz do j = 1, ny do i = 1, nx @@ -147,6 +148,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) press(i,j,k),keyerr,anyerr) else ! use polytropic EOS + keytemp = 0 call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,& rho(i,j,k),eps(i,j,k),xtemp,xye,press(i,j,k),keyerr,anyerr) @@ -189,7 +191,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) z(i,j,k),r(i,j,k),epsnegative,GRHydro_rho_min,pmin, epsmin, & GRHydro_reflevel, GRHydro_C2P_failed(i,j,k), local_perc_ptol) if(GRHydro_C2P_failed(i,j,k).eq.2) then - !$OMP CRITICAL(c2po1) + !$OMP CRITICAL if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then call CCTK_WARN(1,"Convergence problem in c2p") write(warnline,"(A10,i5)") "reflevel: ",GRHydro_reflevel @@ -201,7 +203,7 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) call CCTK_WARN(1,warnline) call CCTK_WARN(0,"Aborting!!!") endif - !$OMP END CRITICAL(c2po1) + !$OMP END CRITICAL endif endif endif @@ -263,6 +265,38 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) #endif end if + + +#if 1 + ! this is needed in the current implementation of refluxing -- in this implementation, + ! the entire correction is applied to the coarse grid cell. + ! This leads to the cell getting too cold. + ! We work around this issue by enforcing that the temperature be + ! above 1.2 MeV at a density above 3.0e11. If this is the case, we inject heat. + if(evolve_temper.ne.0) then + if(GRHydro_eos_hot_eps_fix.ne.0.and.& + rho(i,j,k).gt.5.0d-7 .and. temperature(i,j,k).lt.1.2d0) then + !$OMP CRITICAL +! write(warnline,"(A40)") "Fixing temperature at:" +! call CCTK_WARN(1,warnline) + write(warnline,"(A10,4i6,1P3E15.6)") "fixtemp:", GRHydro_reflevel,i,j,k,& + x(i,j,k),y(i,j,k),z(i,j,k) + call CCTK_WARN(1,warnline) + !$OMP END CRITICAL + keytemp = 1 + temperature(i,j,k) = 1.5d0 + call EOS_Omni_press(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n,& + rho(i,j,k),eps(i,j,k),temperature(i,j,k),y_e(i,j,k),press(i,j,k),keyerr,anyerr) + tau(i,j,k) = sqrt(det) * (rho(i,j,k)+rho(i,j,k)*eps(i,j,k) +press(i,j,k))& + * w_lorentz(i,j,k)**2 - dens(i,j,k) - press(i,j,k) + keytemp = 0 + + endif + endif +#endif + + + end do end do end do @@ -526,7 +560,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve CCTK_REAL epsminl,pminl CCTK_REAL local_perc_ptol - character(len=200) warnline + character(len=256) warnline logical epsnegative integer :: failwarnmode @@ -538,7 +572,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve integer :: nfudgemax,nf n=1;keytemp=0;anyerr=0;keyerr(1)=0 temp0 = 0.0d0;xpress = 0.0d0 - nf=0;nfudgemax=15 + nf=0;nfudgemax=30 ! end EOS Omni vars failinfomode = 1 @@ -553,7 +587,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve y .lt. -1.0d-12 .or.& z .lt. -1.0d-12)) then failwarnmode = 2 - failinfomode = 2 + failinfomode = 2 endif !!$ Undensitize the variables @@ -573,17 +607,8 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve pold = max(1.d-15,xpress) ! error handling if(anyerr.ne.0) then - !$OMP CRITICAL(c2p0) + !$OMP CRITICAL if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then - call CCTK_WARN(failinfomode,"EOS error in c2p 0") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) if(keyerr(1).eq.667.and.GRHydro_eos_hot_eps_fix.ne.0) then ! Handling of the case when no new temperature can be ! found for a given epsilon. The amount of times @@ -592,6 +617,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve nf = 0 do while(anyerr.ne.0.and.nf.le.nfudgemax) anyerr = 0 + utau = utau + rho*abs(epsilon)*0.05d0*w_lorentz**2 epsilon = epsilon + abs(epsilon)*0.05d0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) @@ -602,7 +628,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(failinfomode,warnline) if(nf.gt.nfudgemax) then - call CCTK_WARN(failinfomode,"EOS error in c2p 1: injected heat too many times") + call CCTK_WARN(failinfomode,"EOS error in c2p 0a: injected heat too many times") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(failinfomode,warnline) write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye @@ -614,12 +640,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve call CCTK_WARN(failwarnmode,"Aborting!!!") endif else + call CCTK_WARN(failinfomode,"EOS error in c2p 0") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,dens,epsilon,temp,temp0,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(failinfomode,warnline) call CCTK_WARN(failwarnmode,"Aborting!!!") - endif endif - !$OMP END CRITICAL(c2p0) + !$OMP END CRITICAL endif !!$ Check that the variables have a chance of being physical @@ -652,17 +685,8 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve rho,epsilon,temp,ye,xpress,keyerr,anyerr) ! error handling if(anyerr.ne.0) then - !$OMP CRITICAL(c2p1) + !$OMP CRITICAL if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then - call CCTK_WARN(failinfomode,"EOS error in c2p 1") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) if(keyerr(1).eq.667.and.GRHydro_eos_hot_eps_fix.ne.0) then ! Handling of the case when no new temperature can be ! found for a given epsilon. The amount of times @@ -671,6 +695,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve nf = 0 do while(anyerr.ne.0.and.nf.le.nfudgemax) anyerr = 0 + utau = utau + rho*abs(epsilon)*0.05d0*w_lorentz**2 epsilon = epsilon + abs(epsilon)*0.05d0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) @@ -693,11 +718,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve call CCTK_WARN(failwarnmode,"Aborting!!!") endif else + call CCTK_WARN(failinfomode,"EOS error in c2p 1") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(failinfomode,warnline) call CCTK_WARN(failwarnmode,"Aborting!!!") endif endif - !$OMP END CRITICAL(c2p1) + !$OMP END CRITICAL endif f = pold - xpress @@ -748,25 +781,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve rho,epsilon,temp,ye,xpress,keyerr,anyerr) ! error handling if(anyerr.ne.0) then - !$OMP CRITICAL(c2p2) + !$OMP CRITICAL if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then - call CCTK_WARN(failinfomode,"EOS error in c2p 2") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(failinfomode,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(failinfomode,warnline) if(keyerr(1).eq.667.and.GRHydro_eos_hot_eps_fix.ne.0) then ! Handling of the case when no new temperature can be ! found for a given epsilon. The amount of times ! we get to this place should be monitored, as this ! 'failsafe' mode creates artificial heat. + write(warnline,"(4i5,1P10E15.6)") GRHydro_reflevel,ii,jj,kk,x,y,z + call CCTK_WARN(0,warnline) nf = 0 do while(anyerr.ne.0.and.nf.le.nfudgemax) anyerr = 0 + utau = utau + rho*abs(epsilon)*0.05d0*w_lorentz**2 epsilon = epsilon + abs(epsilon)*0.05d0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& rho,epsilon,temp,ye,xpress,keyerr,anyerr) @@ -789,12 +816,19 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve call CCTK_WARN(failwarnmode,"Aborting!!!") endif else + call CCTK_WARN(failinfomode,"EOS error in c2p 2") + write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(1P10E15.6)") rho,epsilon,temp,ye + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(failinfomode,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(failinfomode,warnline) call CCTK_WARN(failwarnmode,"Aborting!!!") - endif endif - !$OMP END CRITICAL(c2p2) + !$OMP END CRITICAL endif f = pnew - xpress @@ -833,7 +867,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve ! error handling if(anyerr.ne.0) then - !$OMP CRITICAL(c2p3) + !$OMP CRITICAL if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then call CCTK_WARN(failinfomode,"EOS error in c2p 3") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z @@ -846,7 +880,7 @@ subroutine Con2Prim_pt_hot(ii,jj,kk,handle, dens, sx, sy, sz, tau, rho, velx, ve call CCTK_WARN(failinfomode,warnline) call CCTK_WARN(failwarnmode,"Aborting!!!") endif - !$OMP END CRITICAL(c2p3) + !$OMP END CRITICAL endif f = pold - xpress @@ -1767,6 +1801,7 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS) implicit none + DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS integer :: i, j, k, nx, ny, nz @@ -1785,6 +1820,14 @@ subroutine check_GRHydro_C2P_failed(CCTK_ARGUMENTS) if (GRHydro_C2P_failed(i,j,k) == 1) then + if(con2prim_oct_hack.ne.0.and.& + (x(i,j,k) .lt. -1.0d-12 .or.& + y(i,j,k) .lt. -1.0d-12 .or.& + z(i,j,k) .lt. -1.0d-12)) then + cycle + endif + + !$OMP CRITICAL call CCTK_WARN(1,'Con2Prim failed; stopping the code.') call CCTK_WARN(1,'Even with mesh refinement, this point is not restricted from a finer level, so this is really an error') diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 8a428c9..fbe118b 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -254,7 +254,7 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & drho,deps,temp,ye,dpress,keyerr,anyerr) keytemp=0 if(anyerr.ne.0) then - !$OMP CRITICAL(p2c1) + !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(1,warnline) @@ -264,13 +264,13 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & call CCTK_WARN(1,warnline) write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) - !$OMP END CRITICAL(p2c1) + !$OMP END CRITICAL endif else ! This is a way of recovering even on finer refinement levels: ! Use the average temperature at the interface instead of the ! reconstructed specific internal energy. - !$OMP CRITICAL(p2c2) + !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(1,warnline) @@ -280,14 +280,14 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & call CCTK_WARN(1,warnline) write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) - !$OMP END CRITICAL(p2c2) + !$OMP END CRITICAL keytemp=1 temp = temp0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) keytemp=0 if(anyerr.ne.0) then - !$OMP CRITICAL(p2c3) + !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot") write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z call CCTK_WARN(1,warnline) @@ -298,7 +298,7 @@ subroutine prim2con_hot(handle, GRHydro_reflevel, ii, jj, kk, & write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) call CCTK_WARN(0,"Aborting!!!") - !$OMP END CRITICAL(p2c3) + !$OMP END CRITICAL endif endif endif -- cgit v1.2.3