Commits November 2010

commits@lists.einsteintoolkit.org

9 participants
258 discussions

[svn:einsteintoolkit] incoming/EOS_Omni/ (Rev. 23)
by cott＠tapir.caltech.edu 20 Nov '10

20 Nov '10

User: cott Date: 2010/11/20 10:02 PM Modified: /EOS_Omni/ interface.ccl /EOS_Omni/src/ EOS_Omni_MultiVarCalls.F90 /EOS_Omni/src/nuc_eos/ findtemp.F90, nuc_eos.F90 Log: * add a new routine that computes dpde|rho and dpdrho|e at once. File Changes: Directory: /EOS_Omni/ ===================== File [modified]: interface.ccl Delta lines: +16 -0 =================================================================== --- EOS_Omni/interface.ccl 2010-11-17 15:37:53 UTC (rev 22) +++ EOS_Omni/interface.ccl 2010-11-21 04:02:46 UTC (rev 23) @@ -52,7 +52,23 @@ PROVIDES FUNCTION EOS_Omni_DPressByDRho WITH EOS_Omni_EOS_DPressByDRho LANGUAGE Fortran +void FUNCTION EOS_Omni_dpderho_dpdrhoe(CCTK_INT IN eoskey, \ + CCTK_INT IN havetemp, \ + CCTK_REAL IN rf_precision, \ + CCTK_INT IN npoints, \ + CCTK_REAL IN ARRAY rho, \ + CCTK_REAL INOUT ARRAY eps, \ + CCTK_REAL INOUT ARRAY temp, \ + CCTK_REAL IN ARRAY ye, \ + CCTK_REAL OUT ARRAY dpderho, \ + CCTK_REAL OUT ARRAY dpdrhoe, \ + CCTK_INT OUT ARRAY keyerr, \ + CCTK_INT OUT anyerr) +PROVIDES FUNCTION EOS_Omni_dpderho_dpdrhoe WITH EOS_Omni_EOS_dpderho_dpdrhoe LANGUAGE Fortran + + + void FUNCTION EOS_Omni_cs2(CCTK_INT IN eoskey, \ CCTK_INT IN havetemp, \ CCTK_REAL IN rf_precision, \ Directory: /EOS_Omni/src/nuc_eos/ ================================= File [modified]: findtemp.F90 Delta lines: +4 -2 =================================================================== --- EOS_Omni/src/nuc_eos/findtemp.F90 2010-11-17 15:37:53 UTC (rev 22) +++ EOS_Omni/src/nuc_eos/findtemp.F90 2010-11-21 04:02:46 UTC (rev 23) @@ -16,10 +16,10 @@ integer :: rl = 0 integer itmax,i,keyerrt integer ii,jj,kk - + keyerrt=0 - tol=rfeps ! need to find energy to less than 1 in 10^-10 + tol=rfeps ! need to find energy to less than 1 in 10^10 itmax=20 ! use at most 20 iterations, then bomb lt=lt0 @@ -95,12 +95,14 @@ keyerrt=0 goto 12 else +#if 0 ! total failure write(*,*) "EOS: Did not converge in findtemp!" write(*,*) "rl,logrho,logtemp0,ye,lt,eps,eps0,abs(eps-eps0)/eps0" write(*,"(i4,i4,1P10E19.10)") i,rl,lr,lt0,y,lt,eps,eps0,abs(eps-eps0)/eps0 write(*,*) "Tried calling bisection... didn't help... :-/" write(*,*) "Bisection error: ",keyerrt +#endif endif endif File [modified]: nuc_eos.F90 Delta lines: +104 -0 =================================================================== --- EOS_Omni/src/nuc_eos/nuc_eos.F90 2010-11-17 15:37:53 UTC (rev 22) +++ EOS_Omni/src/nuc_eos/nuc_eos.F90 2010-11-21 04:02:46 UTC (rev 23) @@ -348,7 +348,92 @@ end subroutine nuc_eos_press_eps +subroutine nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& + xdpderho,& + keytemp,keyerr,rfeps) + use eosmodule + implicit none + + real*8, intent(in) :: xrho,xye + real*8, intent(inout) :: xtemp,xenr + real*8, intent(in) :: rfeps + real*8, intent(out) :: xdpdrhoe,xdpderho + integer, intent(in) :: keytemp + integer, intent(out) :: keyerr + + ! local variables + real*8 :: xcs2,xprs + real*8 :: lr,lt,y,xx,xeps,leps,xs + real*8 :: d1,d2,d3,ff(8) + integer :: keyerrt = 0 + + if(xrho.gt.eos_rhomax) then + stop "nuc_eos: rho > rhomax" + endif + + if(xrho.lt.eos_rhomin*1.2d0) then + call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + return + endif + + if(xye.gt.eos_yemax) then + stop "nuc_eos: ye > yemax" + endif + + if(xye.lt.eos_yemin) then + keyerr = 102 + write(6,*) "ye: ", xye + stop "nuc_eos ye < yemin" + endif + + if(keytemp.eq.1) then + if(xtemp.gt.eos_tempmax) then + stop "nuc_eos: temp > tempmax" + endif + + if(xtemp.lt.eos_tempmin) then + call nuc_low_eos(xrho,xenr,xprs,xcs2,xdpderho,xdpdrhoe,keytemp) + return + endif + endif + + lr = log10(xrho) + lt = log10(xtemp) + y = xye + xeps = xenr + energy_shift + leps = log10(max(xeps,1.0d0)) + + keyerr = 0 + + if(keytemp.eq.0) then + !need to find temperature based on xeps + call findtemp(lr,lt,y,leps,keyerrt,rfeps) + if(keyerrt.ne.0) then + keyerr = keyerrt + return + endif + xtemp = 10.0d0**lt + + elseif(keytemp.eq.2) then + stop "eos_nuc_press does not support keytemp.eq.2" + endif + + ! have temperature, proceed: + call findall_dpdr_dpde(lr,lt,y,ff) + xdpdrhoe = 10.0d0**ff(1) + xdpderho = 10.0d0**ff(2) + + + if((keytemp.eq.1).or.(keytemp.eq.2)) then + call findthis(lr,lt,y,xeps,alltables(:,:,:,2),d1,d2,d3) + xeps = 10.0d0**xeps - energy_shift + xenr = xeps + endif + +end subroutine nuc_eos_dpdr_dpde + + subroutine findthis(lr,lt,y,value,array,d1,d2,d3) use eosmodule @@ -426,3 +511,22 @@ ff(:) = ffx(:,1) end subroutine findall_press_eps + +subroutine findall_dpdr_dpde(lr,lt,y,ff) + + use eosmodule + implicit none + + real*8 ffx(8,1) + real*8 ff(8) + real*8 lr,lt,y + integer i + integer :: nvarsx = 2 + + +! Ewald's interpolator + call intp3d_many(lr,lt,y,ffx,1,alltables(:,:,:,7:8), & + nrho,ntemp,nye,nvarsx,logrho,logtemp,ye) + ff(:) = ffx(:,1) + +end subroutine findall_dpdr_dpde Directory: /EOS_Omni/src/ ========================= File [modified]: EOS_Omni_MultiVarCalls.F90 Delta lines: +120 -0 =================================================================== --- EOS_Omni/src/EOS_Omni_MultiVarCalls.F90 2010-11-17 15:37:53 UTC (rev 22) +++ EOS_Omni/src/EOS_Omni_MultiVarCalls.F90 2010-11-21 04:02:46 UTC (rev 23) @@ -82,3 +82,123 @@ enddo end subroutine EOS_Omni_EOS_short + + +subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,& + rho,eps,temp,ye,dpderho,dpdrhoe,keyerr,anyerr) + + use EOS_Omni_Module + implicit none + DECLARE_CCTK_PARAMETERS + + CCTK_INT, intent(in) :: eoskey,keytemp,npoints + CCTK_INT, intent(out) :: keyerr(npoints) + CCTK_INT, intent(out) :: anyerr + CCTK_REAL, intent(in) :: rf_precision + CCTK_REAL, intent(in) :: rho(npoints),ye(npoints) + CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints) + CCTK_REAL, intent(out) :: dpderho(npoints) + CCTK_REAL, intent(out) :: dpdrhoe(npoints) + + ! local vars + integer :: i + character(256) :: warnstring + real*8 :: hybrid_local_gamma + real*8 :: hybrid_local_k_cgs + real*8 :: hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2 + ! temporary vars for nuc_eos + real*8 :: xrho,xye,xtemp,xenr,xent + real*8 :: xprs,xmunu,xcs2 + real*8 :: xdedt,xdpderho,xdpdrhoe + + + anyerr = 0 + keyerr(:) = 0 + + select case (eoskey) + case (1) + ! polytropic EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * poly_k_cgs * & + (rho(i)*inv_rho_gf)**(poly_gamma) / & + (poly_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdrhoe(i) = press_gf * poly_k_cgs * & + poly_gamma * inv_rho_gf * & + (rho(i)*inv_rho_gf) ** (poly_gamma - 1.d0) + dpderho(i) = 0.0d0 + enddo + case (2) + ! gamma-law EOS + if(keytemp.eq.1) then + do i=1,npoints + eps(i) = press_gf * gl_k_cgs * & + (rho(i)*inv_rho_gf)**(gl_gamma) / & + (gl_gamma - 1.0d0) / rho(i) + enddo + endif + do i=1,npoints + dpdrhoe(i) = (gl_gamma-1.0d0) * & + eps(i) + dpderho(i) = (gl_gamma - 1.0d0) * & + rho(i) + enddo + case (3) + ! hybrid EOS + do i=1,npoints + if(rho(i).gt.hybrid_rho_nuc) then + hybrid_local_gamma = hybrid_gamma2 + hybrid_local_k_cgs = hybrid_k2_cgs + else + hybrid_local_gamma = hybrid_gamma1 + hybrid_local_k_cgs = hybrid_k1_cgs + endif + hybrid_dp_poly = hybrid_local_gamma * press_gf * & + hybrid_local_k_cgs * rho(i)**(hybrid_local_gamma - 1.0d0) * & + inv_rho_gf**hybrid_local_gamma + + hybrid_dp_th1 = - hybrid_local_gamma * press_gf * hybrid_local_k_cgs * & + (hybrid_gamma_th - 1.d0) / (hybrid_local_gamma - 1.d0) * & + rho(i)**(hybrid_local_gamma - 1.d0) * inv_rho_gf**hybrid_local_gamma + + hybrid_dp_th2 = (hybrid_gamma_th - 1.d0) * eps(i) & + - (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) / & + (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) * & + press_gf * hybrid_k1_cgs * inv_rho_gf**hybrid_gamma1 * & + hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) + + dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2 + dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i) + enddo + case (4) + do i=1,npoints + xrho = rho(i) * inv_rho_gf + xtemp = temp(i) + xye = ye(i) + xenr = eps(i) * inv_eps_gf + call nuc_eos_dpdr_dpde(xrho,xtemp,xye,xenr,xdpdrhoe,& + xdpderho,& + keytemp,keyerr(i),rf_precision) + + if(keyerr(i).ne.0) then + anyerr = 1 + endif + + if(keytemp.eq.1) then + eps(i) = xenr * eps_gf + else + temp(i) = xtemp + endif + + dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf + dpderho(i) = xdpderho * press_gf * inv_eps_gf + enddo + case DEFAULT + write(warnstring,*) "eoskey ",eoskey," not implemented!" + call CCTK_WARN(0,warnstring) + end select + + end subroutine EOS_Omni_EOS_dpderho_dpdrhoe

1 0

[svn:einsteintoolkit] GRHydro/branches/hot_and_MHD_temp_dev/ (Rev. 182)
by cott＠tapir.caltech.edu 20 Nov '10

20 Nov '10

User: cott Date: 2010/11/20 09:46 PM Modified: /branches/hot_and_MHD_temp_dev/ param.ccl, schedule.ccl /branches/hot_and_MHD_temp_dev/src/ GRHydro_Prim2Con.F90, Utils.F90 Log: * modify Prim2Con to handle issues on coarse reflevels when using a hot EOS * update GRHydro_reflevel before FluxTerms File Changes: Directory: /branches/hot_and_MHD_temp_dev/ ========================================== File [modified]: param.ccl Delta lines: +5 -6 =================================================================== --- branches/hot_and_MHD_temp_dev/param.ccl 2010-11-21 03:29:28 UTC (rev 181) +++ branches/hot_and_MHD_temp_dev/param.ccl 2010-11-21 03:46:50 UTC (rev 182) @@ -259,17 +259,11 @@ ####################### Other Parameters ############################## #Parameters controlling conservative <-> primitive change. - int GRHydro_c2p_warnlevel "Warnlevel for Con2Prim warnings" STEERABLE=ALWAYS { 0:1 :: "Either 0 or 1" } 0 -# This parameter is not used anymore in maintained code: the warnings are produced only for really evolved points. -int GRHydro_c2p_warn_from_reflevel "Start warning on given refinement level and on higher levels" STEERABLE=ALWAYS -{ - 0: :: "Greater or equal to 0" -} 0 boolean c2p_reset_pressure "If the pressure guess is unphysical should we recompute it?" { @@ -551,3 +545,8 @@ boolean con2prim_oct_hack "Disregard c2p failures in oct/rotsym90 boundary cells with xyz < 0" STEERABLE=ALWAYS { } "no" + +int GRHydro_c2p_warn_from_reflevel "Start warning on given refinement level and on higher levels" STEERABLE=ALWAYS +{ + 0: :: "Greater or equal to 0" +} 0 File [modified]: schedule.ccl Delta lines: +6 -1 =================================================================== --- branches/hot_and_MHD_temp_dev/schedule.ccl 2010-11-21 03:29:28 UTC (rev 181) +++ branches/hot_and_MHD_temp_dev/schedule.ccl 2010-11-21 03:46:50 UTC (rev 182) @@ -877,8 +877,13 @@ LANG: Fortran } "Calculate current refinement level" +schedule GRHydro_RefinementLevel IN FluxTerms BEFORE Reconstruct +{ + LANG: Fortran +} "Calculate current refinement level" + #debug -###schedule GRHydro_Debug IN ZelmaniLeakChangeYe BEFORE ZelmaniLeak_ye_of_rho +###schedule GRHydro_Debug IN HydroBase_PostStep ###{ ### LANG: Fortran ###} "debug" Directory: /branches/hot_and_MHD_temp_dev/src/ ============================================== File [modified]: GRHydro_Prim2Con.F90 Delta lines: +39 -12 =================================================================== --- branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-21 03:29:28 UTC (rev 181) +++ branches/hot_and_MHD_temp_dev/src/GRHydro_Prim2Con.F90 2010-11-21 03:46:50 UTC (rev 182) @@ -112,7 +112,10 @@ ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the ! variable 'temperature' in prim2con at the interfaces - xtemp = temperature(i,j,k) + ! We will instead use an average temperature as an initial + ! guess. + xtemp = 0.5d0*(temperature(i,j,k) + & + temperature(i-xoffset,j-yoffset,k-zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxl,gxyl,gxzl,gyyl,& gyzl,gzzl, & @@ -126,7 +129,10 @@ ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the ! variable 'temperature' in prim2con at the interfaces - xtemp = temperature(i,j,k) + ! We will instead use an average temperature as an initial + ! guess. + xtemp = 0.5d0*(temperature(i,j,k) + & + temperature(i+xoffset,j+yoffset,k+zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, & i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxr,gxyr,gxzr,& gyyr,gyzr,gzzr, & @@ -238,16 +244,37 @@ drho,deps,temp,ye,dpress,keyerr,anyerr) ! error handling if(anyerr.ne.0) then - 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) - write(warnline,"(1P10E15.6)") drho,deps,temp,ye - call CCTK_WARN(1,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(1,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(1,warnline) - call CCTK_WARN(0,"Aborting!!!") + if(GRHydro_reflevel.lt.GRHydro_c2p_warn_from_reflevel) then + ! in this case (coarse grid error that is hopefully restricted + ! away), we use the average temperature between cells and call + ! the EOS with keytemp=1 + keytemp=1 + 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 + 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) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + endif + else + 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) + write(warnline,"(1P10E15.6)") drho,deps,temp,ye + call CCTK_WARN(1,warnline) + write(warnline,"(A7,i8)") "code: ",keyerr(1) + call CCTK_WARN(1,warnline) + write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel + call CCTK_WARN(1,warnline) + call CCTK_WARN(0,"Aborting!!!") + endif endif vlowx = gxx*dvelx + gxy*dvely + gxz*dvelz File [modified]: Utils.F90 Delta lines: +18 -6 =================================================================== --- branches/hot_and_MHD_temp_dev/src/Utils.F90 2010-11-21 03:29:28 UTC (rev 181) +++ branches/hot_and_MHD_temp_dev/src/Utils.F90 2010-11-21 03:46:50 UTC (rev 182) @@ -32,15 +32,27 @@ implicit none DECLARE_CCTK_ARGUMENTS integer i,j,k + integer nx, ny, nz GRHydro_reflevel = aint(log10(dble(cctk_levfac(1)))/log10(2.0d0)) + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) - i=25 - j=3 - k=3 - write(6,"(i4,1P10E15.6)") GRHydro_reflevel, temperature(i,j,k),temperature(i,j,k),temperature(i,j,k) - write(6,"(i4,1P10E15.6)") GRHydro_reflevel, rho(i,j,k) - write(6,"(3i4)") cctk_lsh(1),cctk_lsh(2),cctk_lsh(3) + if (GRHydro_reflevel .lt. 3) return + + do i=4,nx + do j=4,ny + do k=4,4 + if( r(i,j,k)-4.0d0 .lt. 96.0d0 .and. & + r(i,j,k)+4.0d0 .gt. 96.0d0) then + write(6,"(4i4,1P10E15.6)") GRHydro_reflevel, i,j,k,& + r(i,j,k),rho(i,j,k),eps(i,j,k),y_e(i,j,k),temperature(i,j,k) + endif + enddo + enddo + enddo call flush(6) ! if(GRHydro_reflevel.eq.4.and.temperature(1,1,1).lt.0.1d0) then

1 0