Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 187)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/ param.ccl /trunk/src/ GRHydro_BondiM_new.F90 Log: GRHydro_InitData: Set optionally the plasma beta parameter at the sonic point. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +37 -8 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:26 UTC (rev 186) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:28 UTC (rev 187) @@ -60,12 +60,13 @@ CCTK_REAL :: ONEmTINY, tiny PARAMETER (N_points=100000,ONEmTINY=0.999999d0,tiny=1.0d-12) CCTK_REAL :: M, Msq, Mdot, rs, gam, rmin_bondi, rmax_bondi, cs_sq,cs,vs_sq,vs,rhos,gtemp,hs, Kval, Qdot + CCTK_REAL :: psonic, riso_s CCTK_REAL :: logrmin,dlogr,rhotmp,utmp,vtmp,rspher CCTK_REAL :: r_bondi(N_points), logr_bondi(N_points), rho_bondi(N_points), u_bondi(N_points), v_bondi(N_points) CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso - CCTK_REAL :: bondi_bsmooth + CCTK_REAL :: bondi_bsmooth, bmag character(400) :: debug_message @@ -74,11 +75,17 @@ Msq = M*M Mdot = mdot_sonicpt_bondi rs = r_sonicpt_bondi + riso_s = 0.5d0*(rs-M+sqrt(rs*(rs-2.0d0*M))) gam = gl_gamma - write(debug_message,'(a,4f22.14)') "Bondi_pars: ",M,mdot_sonicpt_bondi, & - r_sonicpt_bondi,gl_gamma + write(debug_message,'(a,2f22.14)') "Bondi_pars: M, mdot_sonicpt_bondi,",& + M,mdot_sonicpt_bondi call CCTK_INFO(debug_message) + write(debug_message,'(a,2f22.14)') "Bondi_pars: r_sonicpt_bondi,gl_gamma", & + r_sonicpt_bondi,gl_gamma + call CCTK_INFO(debug_message) + write(debug_message,'(a,f22.14)') "Bondi_pars: riso_s", riso_s + call CCTK_INFO(debug_message) rmin_bondi = M * bondi_rmin(1) rmax_bondi = M * bondi_rmax(1) @@ -98,13 +105,18 @@ hs = 1. / ( 1. - cs_sq / (gam - 1.) ) Kval = hs * cs_sq * rhos**(-gtemp) / gam Qdot = hs * hs * ( 1. - 3. * vs_sq ) + ! Get the pressure value psonic at the sonic point + psonic = Kval * rhos**gam logrmin = log10(rmin_bondi) dlogr = (log10(rmax_bondi) - logrmin)/(1.*(N_points-1)) - write(debug_message,'(a,8f22.14)') "More pars: ",cs,vs,rhos,hs,Kval,Qdot,& - logrmin,dlogr + write(debug_message,'(a,4f22.14)') "Bondi pars: cs,vs,rhos,hs",cs,vs,rhos,hs call CCTK_INFO(debug_message) + write(debug_message,'(a,2f22.14)') "Bondi pars: Kval,Qdot ",Kval,Qdot + call CCTK_INFO(debug_message) + write(debug_message,'(a,2f22.14)') "Bondi pars: logrmin,dlogr",logrmin,dlogr + call CCTK_INFO(debug_message) rhotmp=1.0d30 imin=1 @@ -122,6 +134,7 @@ !!$ rhotmp = -1. !!$ start with guess rhotmp=rhos !!$ start with value at sonic point! + do i=imin,N_points rspher = r_bondi(i) call find_bondi_solution( rspher, rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot ) @@ -203,6 +216,22 @@ + ! Note that: B^r(t,r) = bondi_bmag M^2 / ( sqrt(\Lambda) \lambda r^2 ) + ! according to eq. 82 of PRD82 084031 (2010). Assuming Schwarzschild + ! coordinates B^r = bondi_bmag M^2/r^2 sqrt(1-2M/r). We can show that + ! b^\mu b_\mu = (B^r)^2 and from the definition of plasma beta parameter + ! we can find that bondi_bmag = sqrt(2P/\beta) r^2/M 1/sqrt(1-2M/r) + if(set_bondi_beta_sonicpt.ne.0) then + bmag = sqrt(2.0d0*psonic/bondi_beta_sonicpt)*rs**2/M & + /sqrt(1.0d0-2.0d0*M/rs) + else + bmag = bondi_bmag + end if + write(debug_message,'(a,2f22.14)')'Bondi pars: rs, psonic',rs, psonic + call CCTK_INFO(debug_message) + write(debug_message,'(a,2f22.14)')'Bondi pars: bondi_bmag,bondi_beta_sonicpt',& + bmag,bondi_beta_sonicpt + call CCTK_INFO(debug_message) do i=1,nx do j=1,ny @@ -261,9 +290,9 @@ det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - Bvecx(i,j,k) = bondi_bsmooth*bondi_bmag*M**2*xhat/sqrt(det)/riso**2 - Bvecy(i,j,k) = bondi_bsmooth*bondi_bmag*M**2*yhat/sqrt(det)/riso**2 - Bvecz(i,j,k) = bondi_bsmooth*bondi_bmag*M**2*zhat/sqrt(det)/riso**2 + Bvecx(i,j,k) = bondi_bsmooth*bmag*M**2*xhat/sqrt(det)/riso**2 + Bvecy(i,j,k) = bondi_bsmooth*bmag*M**2*yhat/sqrt(det)/riso**2 + Bvecz(i,j,k) = bondi_bsmooth*bmag*M**2*zhat/sqrt(det)/riso**2 call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k), & gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +9 -0 =================================================================== --- trunk/param.ccl 2012-11-09 01:59:26 UTC (rev 186) +++ trunk/param.ccl 2012-11-09 01:59:28 UTC (rev 187) @@ -277,6 +277,15 @@ 0:* :: "Anything positive" } 0.01 +CCTK_BOOLEAN set_bondi_beta_sonicpt "Set plasma beta parameter instead of bondi_bmag" +{ +} no + +CCTK_REAL bondi_beta_sonicpt "Plasma beta parameter at the sonic point. Calculate bondi_bmag afterwards." +{ + (0:* :: "positive" +} 1.0 + # For Poloidal Magnetic field test: CCTK_REAL poloidal_A_b "Vector potential strength"

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 186)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ CheckParam.c Log: GRHydro_InitData: fix typo in error message. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: CheckParam.c Delta lines: +1 -1 =================================================================== --- trunk/src/CheckParam.c 2012-11-09 01:59:23 UTC (rev 185) +++ trunk/src/CheckParam.c 2012-11-09 01:59:26 UTC (rev 186) @@ -62,7 +62,7 @@ if((CCTK_Equals(initial_hydro,"magnetized_bondi_solution") || CCTK_Equals(initial_hydro,"magnetized_bondi_solution_iso")) && !CCTK_Equals(initial_Bvec, "magnetized Bondi") ) { - CCTK_PARAMWARN("Please set initial_Bec='magnetized Bondi' to properly initialize the magnetic field for Bondi solutions!"); + CCTK_PARAMWARN("Please set initial_Bvec='magnetized Bondi' to properly initialize the magnetic field for Bondi solutions!"); } if( CCTK_Equals(initial_hydro,"magnetized_bondi_solution") || CCTK_Equals(initial_hydro,"hydro_bondi_solution") ) {

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 185)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/ interface.ccl /trunk/src/ GRHydro_C2P2CM.F90, GRHydro_P2C2PM.F90 Log: ng the interface for magentized con2prim for consistency From: Josh Faber <jfaberrit(a)jfaber3.local> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_C2P2CM.F90 Delta lines: +4 -3 =================================================================== --- trunk/src/GRHydro_C2P2CM.F90 2012-11-09 01:59:21 UTC (rev 184) +++ trunk/src/GRHydro_C2P2CM.F90 2012-11-09 01:59:23 UTC (rev 185) @@ -90,7 +90,7 @@ sy_send = 0.166666658d0 sz_send = 0.166666658d0 tau_send = 0.484123939d0 - + bvcx_send = Bx_init bvcy_send = By_init bvcz_send = Bz_init @@ -98,6 +98,7 @@ bconsy_send = By_init bconsz_send = Bz_init + rho_send=1.0d-6 eps_send = 1.0d-6 press_send = 6.666666666666667d-7 w_lorentz_send = 1.0d0 @@ -138,8 +139,8 @@ write(*,*) ' Bvecz : ',bvcz_send write(*,*) 'C2P2CM test: getting the associated primitive variables.' - call Con2PrimGenM(GRHydro_eos_handle,local_gam(1),dens_send,sx_send,sy_send,sz_send, & - tau_send,bconsx_send,bconsy_send,bconsz_send,rho_send,velx_send,vely_send,velz_send, & + call Con2PrimGenM(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & + tau_send,bconsx_send,bconsy_send,bconsz_send,xtemp(1),xye(1),rho_send,velx_send,vely_send,velz_send, & eps_send,press_send,bvcx_send,bvcy_send,bvcz_send,b2_send,w_lorentz_send, & gxx_send,gxy_send,gxz_send,gyy_send,gyz_send,gzz_send,& uxx,uxy,uxz,uyy,uyz,uzz,det,& File [modified]: GRHydro_P2C2PM.F90 Delta lines: +2 -2 =================================================================== --- trunk/src/GRHydro_P2C2PM.F90 2012-11-09 01:59:21 UTC (rev 184) +++ trunk/src/GRHydro_P2C2PM.F90 2012-11-09 01:59:23 UTC (rev 185) @@ -153,9 +153,9 @@ write(*,*) 'P2C2PM test: getting the associated primitive variables.' - call Con2PrimGenM(GRHydro_eos_handle,local_gam(1),dens_send,sx_send,sy_send,sz_send, & + call Con2PrimGenM(GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,local_gam(1),dens_send,sx_send,sy_send,sz_send, & tau_send,bconsx_send,bconsy_send,bconsz_send,& - rho_send(1),velx_send,vely_send,velz_send, & + xtemp(1),xye(1),rho_send(1),velx_send,vely_send,velz_send, & eps_send(1),press_send(1), & bvcx_send,bvcy_send,bvcz_send,b2_send,& w_lorentz_send, & Directory: /trunk/ ================== File [modified]: interface.ccl Delta lines: +2 -19 =================================================================== --- trunk/interface.ccl 2012-11-09 01:59:21 UTC (rev 184) +++ trunk/interface.ccl 2012-11-09 01:59:23 UTC (rev 185) @@ -114,10 +114,10 @@ CCTK_INT INOUT GRHydro_reflevel, \ CCTK_REAL OUT retval) -void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ +void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ - CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ @@ -131,23 +131,6 @@ CCTK_INT OUT epsnegative, \ CCTK_REAL OUT retval) -#void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_INT INOUT keytemp, CCTK_REAL INOUT prec, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ -# CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ -# CCTK_REAL INOUT tau, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ -# CCTK_REAL INOUT y_e, CCTK_REAL INOUT temp, CCTK_REAL INOUT rho, \ -# CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ -# CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ -# CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ -# CCTK_REAL INOUT Bvecsq, \ -# CCTK_REAL INOUT w_lorentz, \ -# CCTK_REAL INOUT gxx, CCTK_REAL INOUT gxy, CCTK_REAL INOUT gxz, \ -# CCTK_REAL INOUT gyy, CCTK_REAL INOUT gyz, CCTK_REAL INOUT gzz, \ -# CCTK_REAL INOUT uxx, CCTK_REAL INOUT uxy, CCTK_REAL INOUT uxz, \ -# CCTK_REAL INOUT uyy, CCTK_REAL INOUT uyz, CCTK_REAL INOUT uzz, \ -# CCTK_REAL INOUT det, \ -# CCTK_INT OUT epsnegative, \ -# CCTK_REAL OUT retval) - void FUNCTION Con2PrimPolyM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 184)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_BondiM_new.F90, GRHydro_Bondi_new.F90 Log: GRHydro_Init_Data: only write bondi data file on root process From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +8 -6 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:19 UTC (rev 183) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:21 UTC (rev 184) @@ -149,12 +149,14 @@ v_bondi(i) = vtmp enddo - open (47,file="bondi.asc",form="formatted") - do i=1,N_points - write(47,'(i5,4f22.14)')i,r_bondi(i),rho_bondi(i),& - u_bondi(i),v_bondi(i) - end do - close(47) + if (CCTK_MyProc(cctkGH) == 0) then + open (47,file="bondi.asc",form="formatted") + do i=1,N_points + write(47,'(i5,4f22.14)')i,r_bondi(i),rho_bondi(i),& + u_bondi(i),v_bondi(i) + end do + close(47) + end if !!$ write(debug_message,'(a,4f22.14)') "i=1:",r_bondi(1),rho_bondi(1),& !!$ u_bondi(1),v_bondi(1) File [modified]: GRHydro_Bondi_new.F90 Delta lines: +8 -6 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:19 UTC (rev 183) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:21 UTC (rev 184) @@ -142,12 +142,14 @@ v_bondi(i) = vtmp enddo - open (47,file="bondi.asc",form="formatted") - do i=1,N_points - write(47,'(i5,4f22.14)')i,r_bondi(i),rho_bondi(i),& - u_bondi(i),v_bondi(i) - end do - close(47) + if(CCTK_MyProc(cctkGH) .eq. 0) then + open (47,file="bondi.asc",form="formatted") + do i=1,N_points + write(47,'(i5,4f22.14)')i,r_bondi(i),rho_bondi(i),& + u_bondi(i),v_bondi(i) + end do + close(47) + end if !!$ write(debug_message,'(a,4f22.14)') "i=1:",r_bondi(1),rho_bondi(1),& !!$ u_bondi(1),v_bondi(1)

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 183)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/ param.ccl /trunk/src/ CheckParam.c Log: GRHydro_Init_Data: add "magnetized Bondi" to hydrobase::initial_bvec this is so that there is a proper initial_bvec setting rather than just "zero". From: Roland Haas <roland.haas(a)physics.gatech.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: CheckParam.c Delta lines: +8 -2 =================================================================== --- trunk/src/CheckParam.c 2012-11-09 01:59:16 UTC (rev 182) +++ trunk/src/CheckParam.c 2012-11-09 01:59:19 UTC (rev 183) @@ -50,15 +50,21 @@ if(CCTK_Equals(Bvec_evolution_method,"GRHydro") && CCTK_Equals(initial_hydro,"hydro_bondi_solution") ) { - CCTK_PARAMWARN("Please use set initial_hydro=magnetized_bondi_solution instead to initialize the magnetic field correctly!"); + CCTK_PARAMWARN("Please set initial_hydro='magnetized_bondi_solution' instead to initialize the magnetic field correctly!"); } if(!CCTK_Equals(Bvec_evolution_method,"GRHydro") && CCTK_Equals(initial_hydro,"magnetized_bondi_solution") ) { - CCTK_PARAMWARN("Please use set initial_hydro=hydro_bondi_solution instead to NOT initialize the magnetic field!"); + CCTK_PARAMWARN("Please set initial_hydro='hydro_bondi_solution' instead to NOT initialize the magnetic field!"); } + if((CCTK_Equals(initial_hydro,"magnetized_bondi_solution") || CCTK_Equals(initial_hydro,"magnetized_bondi_solution_iso")) && + !CCTK_Equals(initial_Bvec, "magnetized Bondi") ) + { + CCTK_PARAMWARN("Please set initial_Bec='magnetized Bondi' to properly initialize the magnetic field for Bondi solutions!"); + } + if( CCTK_Equals(initial_hydro,"magnetized_bondi_solution") || CCTK_Equals(initial_hydro,"hydro_bondi_solution") ) { if( num_bondi_sols > 1 ) { CCTK_PARAMWARN("Currently only one Bondi solution is supported, please change [num_bondi_sols] "); Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +1 -0 =================================================================== --- trunk/param.ccl 2012-11-09 01:59:16 UTC (rev 182) +++ trunk/param.ccl 2012-11-09 01:59:19 UTC (rev 183) @@ -31,6 +31,7 @@ "shocktube" :: "Shocktube type" "cylexp" :: "Poloidal Magnetic Field" "poloidalmagfield" :: "Poloidal Magnetic Field" + "magnetized Bondi" :: "radial magnetic field appropriate for Bondi test" } shares:ADMBase

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 182)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_BondiM_new.F90 Log: GRHydro_InitData: yet another initialization bug fixed in bondiM. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +2 -1 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:14 UTC (rev 181) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:16 UTC (rev 182) @@ -201,7 +201,6 @@ - bondi_bsmooth = 1.0d0 do i=1,nx do j=1,ny @@ -216,6 +215,8 @@ yhat = yp/riso zhat = zp/riso roverm = riso/M + + bondi_bsmooth = 1.0d0 if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 181)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_BondiM_new.F90 Log: GRHydro_InitData: Fix B-field smoothing in BondiM routine. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +9 -10 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:11 UTC (rev 180) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:14 UTC (rev 181) @@ -65,6 +65,7 @@ CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso + CCTK_REAL :: bondi_bsmooth character(400) :: debug_message @@ -199,11 +200,9 @@ nz = cctk_lsh(3) - write(debug_message,'(a,3f22.14)') 'Lower bound coordinates:',x(1,1,1),y(1,1,1),z(1,1,1) - call CCTK_INFO(debug_message) - write(debug_message,'(a,3f22.14)') 'Upper bound coordinates:',x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz) - call CCTK_INFO(debug_message) - + + bondi_bsmooth = 1.0d0 + do i=1,nx do j=1,ny do k=1,nz @@ -238,7 +237,7 @@ rho(i,j,k) = rhocheck+drhodr*riso*(riso-M)/M else rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 - bondi_bmag = bondi_bmag*8.0d0*riso**3 + bondi_bsmooth = 8.0d0*riso**3 endif utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.) @@ -259,10 +258,10 @@ det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - Bvecx(i,j,k) = bondi_bmag*M**2*xhat/sqrt(det)/riso**2 - Bvecy(i,j,k) = bondi_bmag*M**2*yhat/sqrt(det)/riso**2 - Bvecz(i,j,k) = bondi_bmag*M**2*zhat/sqrt(det)/riso**2 - + Bvecx(i,j,k) = bondi_bsmooth*bondi_bmag*M**2*xhat/sqrt(det)/riso**2 + Bvecy(i,j,k) = bondi_bsmooth*bondi_bmag*M**2*yhat/sqrt(det)/riso**2 + Bvecz(i,j,k) = bondi_bsmooth*bondi_bmag*M**2*zhat/sqrt(det)/riso**2 + call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k), & gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), & det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k), &

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 180)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_BondiM_new.F90 Log: GRHydro_InitData: Force mag field to go to zero at the puncture. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +1 -0 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:09 UTC (rev 179) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:11 UTC (rev 180) @@ -238,6 +238,7 @@ rho(i,j,k) = rhocheck+drhodr*riso*(riso-M)/M else rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 + bondi_bmag = bondi_bmag*8.0d0*riso**3 endif utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.)

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 179)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_BondiM_new.F90, GRHydro_Bondi_new.F90 Log: GRHydro_InitData: port smoother velocity matching to mag. bondi routine. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_BondiM_new.F90 Delta lines: +41 -12 =================================================================== --- trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:06 UTC (rev 178) +++ trunk/src/GRHydro_BondiM_new.F90 2012-11-09 01:59:09 UTC (rev 179) @@ -64,7 +64,7 @@ CCTK_REAL :: r_bondi(N_points), logr_bondi(N_points), rho_bondi(N_points), u_bondi(N_points), v_bondi(N_points) CCTK_REAL :: drhodr, det, rhocheck, rhocheck2, riso, rnew, rsch, ucheck CCTK_REAL :: uiso, uisocheck, vcheck, ucheck2, vcheck2, xhat,yhat, zhat, xp, yp, zp - CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm + CCTK_REAL :: f,df,ddf,a,b,c,rsm,roverm,dudr,uisocheck2,auiso,buiso character(400) :: debug_message @@ -118,7 +118,8 @@ endif enddo - rhotmp = -1. !!$ start with guess +!!$ rhotmp = -1. !!$ start with guess + rhotmp=rhos !!$ start with value at sonic point! do i=imin,N_points rspher = r_bondi(i) @@ -132,7 +133,8 @@ v_bondi(i) = vtmp end do - rhotmp = -1. +!!$ rhotmp = -1. + rhotmp=rhos !!$ start with value at sonic point! do i=imin-1,1,-1 rspher = r_bondi(i) @@ -165,27 +167,42 @@ !!$ write(debug_message,'(a,4f22.14)') "i=1500:",r_bondi(1500),rho_bondi(1500),& !!$ u_bondi(1500),v_bondi(1500) !!$ call CCTK_INFO(debug_message) - -!!$ // find the derivative near r=M +!!$ // find the derivative near r=M in isotropic coords = r=9/4M in schwarzschild; rnew = 2.25 * M - j = floor ( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor ((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck = rho_bondi(j) call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot ) uisocheck = 4.0*vcheck/3.0 +!!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic rnew = 0.25 * 3.02**2 * M/1.01 - j = floor( 0.5 + (log10(rnew) - logrmin) / dlogr ) + j = floor((log10(rnew) - logrmin) / dlogr + 1.0) +!!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck2 = rho_bondi(j) call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot ) + uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02) + drhodr = 100.0*(rhocheck2-rhocheck)/M +!!$ Don't divide by M here, to simplify the math + dudr = 100.0*(uisocheck2-uisocheck) + write(debug_message,'(a,3f22.14)') 'Rhocheck:',rhocheck,rhocheck2,drhodr call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Ucheck:',uisocheck,uisocheck2,dudr + call CCTK_INFO(debug_message) nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) + + + write(debug_message,'(a,3f22.14)') 'Lower bound coordinates:',x(1,1,1),y(1,1,1),z(1,1,1) + call CCTK_INFO(debug_message) + write(debug_message,'(a,3f22.14)') 'Upper bound coordinates:',x(nx,ny,nz),y(nx,ny,nz),z(nx,ny,nz) + call CCTK_INFO(debug_message) do i=1,nx do j=1,ny @@ -203,12 +220,17 @@ if(roverm > ONEmTINY) then rsch = 0.25 * ( 2.*riso + M)**2 / riso - jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) +!!$ jb = floor( 0.5 + (log10(rsch) - logrmin) / dlogr ) + jb = floor( (log10(rsch) - logrmin) / dlogr + 1.0) +!!$ jb = NINT((log10(rsch) - logrmin) / dlogr + 1.0) - if(jb > N_points)jb = N_points + if(jb >= N_points) jb = N_points-1 - rhotmp = rho_bondi(jb) - call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) + rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*& + (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) +!!$ call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(jb),rho_bondi(jb+1)) + call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else @@ -218,7 +240,14 @@ rho(i,j,k) = (rhocheck-drhodr*M/4.0)*(1.-cos(2.*M_PI*riso/M))/2.0 endif utmp = Kval * rho(i,j,k)**( gam ) / (gam - 1.) - uiso = uisocheck * riso / M + +!!$ match to uiso and dudr at roverm=1 +!!$ a R + b R^3 ---> a+b = uisocheck; a+3b = dudr +!!$ b = (dudr-uisocheck)/2; a=3*uisocheck-dudr)/2 + + auiso = 1.5*uisocheck-0.5*dudr + buiso = 0.5*dudr-0.5*uisocheck + uiso = roverm*(auiso+buiso*roverm**2) endif eps(i,j,k) = utmp/rhotmp File [modified]: GRHydro_Bondi_new.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:06 UTC (rev 178) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:09 UTC (rev 179) @@ -148,7 +148,7 @@ u_bondi(i),v_bondi(i) end do close(47) - + !!$ write(debug_message,'(a,4f22.14)') "i=1:",r_bondi(1),rho_bondi(1),& !!$ u_bondi(1),v_bondi(1) !!$ call CCTK_INFO(debug_message)

1 0

[svn:einsteintoolkit] GRHydro_InitData/trunk/src/ (Rev. 178)
by rhaas＠tapir.caltech.edu 08 Nov '12

08 Nov '12

User: rhaas Date: 2012/11/08 07:59 PM Modified: /trunk/src/ GRHydro_Bondi_new.F90 Log: GRHydro_InitData: exchange find_bondi_solution_bracket for find_bondi_solution and improve smoothing of v field. From: Bruno Coutinho Mundim <bcmsma(a)astro.rit.edu> File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Bondi_new.F90 Delta lines: +10 -7 =================================================================== --- trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:02 UTC (rev 177) +++ trunk/src/GRHydro_Bondi_new.F90 2012-11-09 01:59:06 UTC (rev 178) @@ -162,13 +162,14 @@ !!$ u_bondi(1500),v_bondi(1500) !!$ call CCTK_INFO(debug_message) -!!$ // find the derivative near r=M +!!$ // find the derivative near r=M in isotropic coords = r=9/4M in schwarzschild; rnew = 2.25 * M j = floor ((log10(rnew) - logrmin) / dlogr + 1.0) !!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck = rho_bondi(j) - call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & - rho_bondi(j),rho_bondi(j+1)) +!!$ call find_bondi_solution_bracket(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(j),rho_bondi(j+1)) + call find_bondi_solution(rnew,rhocheck, ucheck, vcheck, rs, rhos, M, Mdot, Kval, gam, Qdot) uisocheck = 4.0*vcheck/3.0 !!$ the previous point was r=M in isotropic coords = r=9/4M in schwarzschild; this one is r=1.01M in isotropic @@ -176,8 +177,9 @@ j = floor((log10(rnew) - logrmin) / dlogr + 1.0) !!$ j = NINT((log10(rnew) - logrmin) / dlogr + 1.0) rhocheck2 = rho_bondi(j) - call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & - rho_bondi(j),rho_bondi(j+1)) +!!$ call find_bondi_solution_bracket( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(j),rho_bondi(j+1)) + call find_bondi_solution( rnew, rhocheck2, ucheck2, vcheck2, rs, rhos, M, Mdot, Kval, gam, Qdot) uisocheck2 = vcheck2 / (1.0 - 1.0/2.02) / (1.0+ 1.0/2.02) drhodr = 100.0*(rhocheck2-rhocheck)/M @@ -224,8 +226,9 @@ rhotmp = rho_bondi(jb)+(rho_bondi(jb+1)-rho_bondi(jb))*& (rsch-r_bondi(jb))/(r_bondi(jb+1)-r_bondi(jb)) - call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & - rho_bondi(jb),rho_bondi(jb+1)) +!!$ call find_bondi_solution_bracket( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot, & +!!$ rho_bondi(jb),rho_bondi(jb+1)) + call find_bondi_solution( rsch,rhotmp, utmp, vtmp, rs, rhos, M, Mdot, Kval, gam, Qdot) rho(i,j,k) = rhotmp uiso = vtmp / (1.0 - M/2.0/riso) / (1.0+ M/2.0/riso) else

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits