Commits November 2012

commits@lists.einsteintoolkit.org

9 participants
90 discussions

[svn:einsteintoolkit] WeylScal4/trunk/src/ (Rev. 124)
by schnetter＠cct.lsu.edu 13 Nov '12

13 Nov '12

User: eschnett Date: 2012/11/13 11:44 AM Modified: /trunk/src/ Differencing.h, invars_calc_2nd.cc, invars_calc_4th.cc, invars_calc_Nth.cc, psi4_calc_2nd.cc, psi4_calc_4th.cc, psi4_calc_Nth.cc, psis_calc_2nd.cc, psis_calc_4th.cc, psis_calc_Nth.cc Log: Regenerate code File Changes: Directory: /trunk/src/ ====================== File [modified]: Differencing.h Delta lines: <Unavailable> File [modified]: invars_calc_2nd.cc Delta lines: <Unavailable> File [modified]: invars_calc_4th.cc Delta lines: <Unavailable> File [modified]: invars_calc_Nth.cc Delta lines: <Unavailable> File [modified]: psi4_calc_2nd.cc Delta lines: <Unavailable> File [modified]: psi4_calc_4th.cc Delta lines: <Unavailable> File [modified]: psi4_calc_Nth.cc Delta lines: <Unavailable> File [modified]: psis_calc_2nd.cc Delta lines: <Unavailable> File [modified]: psis_calc_4th.cc Delta lines: <Unavailable> File [modified]: psis_calc_Nth.cc Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] pyGWAnalysis/trunk/ (Rev. 15)
by reisswig＠tapir.caltech.edu 12 Nov '12

12 Nov '12

User: reisswig Date: 2012/11/12 05:00 AM Modified: /trunk/DataAnalysis/ DataAnalysis.py, RadiatedQuantities.py /trunk/DiscreteFunction/ DiscreteFunction.py, WaveFunction.py Log: 1) Radiated energy flux from RWZM. 2) Some operators for DiscreteFunction 3) Load mode array by specifying file format. File Changes: Directory: /trunk/DataAnalysis/ =============================== File [modified]: DataAnalysis.py Delta lines: +21 -5 =================================================================== --- trunk/DataAnalysis/DataAnalysis.py 2012-05-11 00:23:43 UTC (rev 14) +++ trunk/DataAnalysis/DataAnalysis.py 2012-11-12 11:00:27 UTC (rev 15) @@ -427,20 +427,36 @@ res[-1].append([]) return res -def LoadModeArray(basename, lmax, col_time=1, col_re=2, col_im=3): +def LoadModeArray(basename, lmax, col_time=1, col_re=2, col_im=3, filenameformat="cce_conversionscript"): """Loads an array of modes up to lmax from disk. We assume that files begin with 'basename' and end with _l*m*.dat""" res = InitModeArray(lmax) for ll in range(2, lmax+1): for mm in range(-ll, ll+1): - res[ll][mm] = WaveFunction([], []) - res[ll][mm].Load(basename+"_l"+str(ll)+"m"+str(mm)+".dat", col_time, col_re, col_im) - + if (filenameformat == "cce_conversion_script"): + res[ll][mm] = WaveFunction([], []) + res[ll][mm].Load(basename+"_l"+str(ll)+"m"+str(mm)+".dat", col_time, col_re, col_im) + continue + if (filenameformat == "pittnull"): + res[ll][-mm] = WaveFunction([], []) + # Get spaces from Fortran outpout filenames right + tmp1 = "%02d" % ll + if mm < 0: + tmp2 = "m" + else: + tmp2 = "p" + tmp2 += "%02d" % abs(mm) + # Take into account that nullcode has a weird convention for psi4! + res[ll][-mm].Load(basename+".L"+tmp1+"M"+tmp2+".asc", col_time, col_re, col_im) + res[ll][-mm].f *= 2.0*(-1.0)**(mm+1) + res[ll][-mm].f = map(lambda x: x.conjugate(), res[ll][-mm].f) + continue + print("Error loading mode array: filenameformat not recognized!") + return res - def ReconFromsYlm(arrayWF, theta, phi, lmax): """Given an array of WF modes h_lm reconstruct h at angle theta, phi up to l=lmax""" res = WaveFunction(arrayWF[2][0].x, arrayWF[2][0].f*0.0) File [modified]: RadiatedQuantities.py Delta lines: +24 -0 =================================================================== --- trunk/DataAnalysis/RadiatedQuantities.py 2012-05-11 00:23:43 UTC (rev 14) +++ trunk/DataAnalysis/RadiatedQuantities.py 2012-11-12 11:00:27 UTC (rev 15) @@ -39,6 +39,30 @@ +def RadiatedEnergyFluxRWZM(Qeven, Qodd, lmax): + """ Given the time integral of rPsi4 (l,m) modes (a DiscreteFunction mode array), + this function computes the radiated energy flux + that is emitted as a function of time in units of [M].""" + + npoints = Qeven[2][0].Length() + Edot = DiscreteFunction(Qeven[2][0].x, zeros(npoints, float)) + + dtQeven = InitModeArray(2) + for l in range(2, lmax+1): + for m in range(-l, l+1): + dtQeven[l][m] = Qeven[l][m].DifferentiateFunction() + + for ii in range(Edot.Length()): + e = 0.0 + for l in range(2, lmax+1): + for m in range(-l, l+1): + e = e + abs(dtQeven[l][m].f[ii])**2 + abs(Qodd[l][m].f[ii])**2 + Edot.f[ii] = e/(32.*pi) + + return Edot + + + def RadiatedEnergyFlux(int_rPsi4_ModeArray, lmax): """ Given the time integral of rPsi4 (l,m) modes (a DiscreteFunction mode array), this function computes the radiated energy flux Directory: /trunk/DiscreteFunction/ =================================== File [modified]: DiscreteFunction.py Delta lines: +68 -2 =================================================================== --- trunk/DiscreteFunction/DiscreteFunction.py 2012-05-11 00:23:43 UTC (rev 14) +++ trunk/DiscreteFunction/DiscreteFunction.py 2012-11-12 11:00:27 UTC (rev 15) @@ -32,10 +32,13 @@ return False - def __init__(self, x_, f_): + def __init__(self, x_, f_, g_=None): """Initialise a 1D transformation function""" self.x = array(x_) - self.f = array(f_) + if g_ != None: + self.f = array(f_) + 1j*array(g_) + else: + self.f = array(f_) if self.HasUniformSpacing() == True: self.dx = self.x[1]-self.x[0] else: @@ -383,12 +386,31 @@ return xroot + def DifferentiateFunction(self): + """ Differentiate this function and return a function""" + + npoints = self.Length() + res = self.__class__(self.x, self.f*0.0) + + # TODO: fail if not uniform spacing + self.dx = self.x[1]-self.x[0] + res.dx = self.dx + + # Simple second order centred differentiation + + for ii in range (1, self.Length()-1): + res.f[ii] = (self.f[ii+1] - self.f[ii-1])/(2*res.dx) + + return res + + def IntegrateFunction(self): """ Integrate this function and return a function""" npoints = self.Length() res = self.__class__(self.x, self.f*0.0) + # TODO: fail if not uniform spacing self.dx = self.x[1]-self.x[0] res.dx = self.dx @@ -572,6 +594,12 @@ return res + def __add__(self, rval): + res = self.Union(rval) + for ii in range(0, res.Length()): + res.f[ii] = self.Interpolate(res.x[ii]) + rval.Interpolate(res.x[ii]) + return res + def __sub__(self, rval): res = self.Union(rval) for ii in range(0, res.Length()): @@ -591,7 +619,45 @@ res.f[ii] = self.Interpolate(res.x[ii]) / rval.Interpolate(res.x[ii]) return res + def __abs__(self): + npoints = len(self.x) + res = self.__class__(self.x, zeros(npoints, float)) + for ii in range(0, Length()): + res.f[ii] = abs(self.f[ii]) + return res + def __pos__(self): + npoints = len(self.x) + res = self.__class__(self.x, self.f) + return res + + def __neg__(self): + npoints = len(self.x) + res = self.__class__(self.x, -self.f) + return res + + def real(self): + npoints = len(self.x) + res = self.__class__(self.x, zeros(npoints, float)) + for ii in range(0, self.Length()): + res.f[ii] = self.f[ii].real + return res + + def imag(self): + npoints = len(self.x) + res = self.__class__(self.x, zeros(npoints, float)) + for ii in range(0, self.Length()): + res.f[ii] = self.f[ii].imag + return res + + def conjugate(self): + npoints = len(self.x) + res = self.__class__(self.x, zeros(npoints, type(self.f[0]))) + for ii in range(0, self.Length()): + res.f[ii] = self.f[ii].conjugate() + return res + + # implement iterator def __iter__(self): self.index = 0 File [modified]: WaveFunction.py Delta lines: +5 -2 =================================================================== --- trunk/DiscreteFunction/WaveFunction.py 2012-05-11 00:23:43 UTC (rev 14) +++ trunk/DiscreteFunction/WaveFunction.py 2012-11-12 11:00:27 UTC (rev 15) @@ -19,10 +19,13 @@ class WaveFunction(DiscreteFunction): """Class describing discrete 1d wave functions""" - def __init__(self, x_, f_): + def __init__(self, x_, f_, g_=None): """Initialise a 1D wave function""" self.x = array(x_) - self.f = array(f_) + if g_ != None: + self.f = array(f_) + 1j*array(g_) + else: + self.f = array(f_) if self.HasUniformSpacing() == True: self.dx = self.x[1]-self.x[0] else:

1 0

[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

← Newer
1
2
3
4
5
6
7
8
9
Older →

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits November 2012