Commits April 2011

commits@lists.einsteintoolkit.org

12 participants
208 discussions

[svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 123)
by bcmsma＠astro.rit.edu 27 Apr '11

27 Apr '11

User: bmundim Date: 2011/04/27 05:25 PM Modified: /trunk/ param.ccl, schedule.ccl /trunk/par/ balsara1a.par /trunk/src/ CheckParam.c, GRHydro_ShockTube.F90, GRHydro_ShockTubeM.F90, make.code.defn Log: RIT MHD dev: Balsara's and Komissarov's shock tube tests. File Changes: Directory: /trunk/par/ ====================== File [modified]: balsara1a.par Delta lines: +12 -6 =================================================================== --- trunk/par/balsara1a.par 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/par/balsara1a.par 2011-04-27 22:25:23 UTC (rev 123) @@ -5,7 +5,8 @@ driver::ghost_size=3 grhydro::grhydro_stencil=3 -time::dtfac = 0.25 +#time::dtfac = 0.25 +time::dtfac = 0.8 methodoflines::ODE_Method = "rk2" methodoflines::MoL_Intermediate_Steps=2 @@ -37,25 +38,30 @@ grid::type = "BySpacing" grid::domain = "full" -grid::dxyz = 0.01 +#grid::dxyz = 0.01 +grid::dxyz = 0.000625 -driver::global_nx = 150 +driver::global_nx = 1600 driver::global_ny = 10 driver::global_nz = 10 +Cactus::terminate="time" +Cactus::cctk_final_time = 0.4 #cactus::cctk_itlast = 40 -cactus::cctk_itlast = 200 +#cactus::cctk_itlast = 200 +#cactus::cctk_itlast = 2 -IO::out_dir = "balsara1a" +IO::out_dir = $parfile #IOBasic::outInfo_every = 1 #IOBasic::outInfo_vars = "HydroBase::rho" CarpetIOBasic::outInfo_vars="hydrobase::rho" CarpetIOBasic::outInfo_every=1 CarpetIOASCII::out1D_criterion = "time" CarpetIOASCII::out1D_dt = 0.01 +CarpetIOASCII::out1D_d=no CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" CarpetIOHDF5::out_every = 1 -CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" +CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel HydroBase::w_lorentz grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" Directory: /trunk/src/ ====================== File [modified]: CheckParam.c Delta lines: +18 -7 =================================================================== --- trunk/src/CheckParam.c 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/CheckParam.c 2011-04-27 22:25:23 UTC (rev 123) @@ -18,20 +18,31 @@ (CCTK_Equals(initial_hydro,"simple_wave")) || (CCTK_Equals(initial_data,"con2primtest")) || (CCTK_Equals(initial_data,"reconstruction_test")) || - (CCTK_Equals(shocktube_type,"diagshock")) || (CCTK_Equals(shocktube_type,"sphere")))) { CCTK_PARAMWARN("That test not yet implemented in MHD!"); } if(!CCTK_Equals(Bvec_evolution_method,"GRHydro") && - ((CCTK_Equals(shock_case,"Balsara1")) || - (CCTK_Equals(shock_case,"Balsara2")) || - (CCTK_Equals(shock_case,"Balsara3")) || - (CCTK_Equals(shock_case,"Balsara4")) || - (CCTK_Equals(shock_case,"Balsara5")))) + ((CCTK_Equals(shock_case,"Balsara0" )) || + (CCTK_Equals(shock_case,"Balsara1" )) || + (CCTK_Equals(shock_case,"Balsara2" )) || + (CCTK_Equals(shock_case,"Balsara3" )) || + (CCTK_Equals(shock_case,"Balsara4" )) || + (CCTK_Equals(shock_case,"Balsara5" )) || + (CCTK_Equals(shock_case,"Alfven" )) || + (CCTK_Equals(shock_case,"Komissarov1")) || + (CCTK_Equals(shock_case,"Komissarov2")) || + (CCTK_Equals(shock_case,"Komissarov3")) || + (CCTK_Equals(shock_case,"Komissarov4")) || + (CCTK_Equals(shock_case,"Komissarov5")) || + (CCTK_Equals(shock_case,"Komissarov6")) || + (CCTK_Equals(shock_case,"Komissarov7")) || + (CCTK_Equals(shock_case,"Komissarov8")) || + (CCTK_Equals(shock_case,"Komissarov9")) || + (CCTK_Equals(initial_hydro,"cylexp")) + )) { CCTK_PARAMWARN("That test requires MHD! Set Bvec_evolution_method=GRHYDRO!"); } } - File [modified]: GRHydro_ShockTube.F90 Delta lines: +11 -0 =================================================================== --- trunk/src/GRHydro_ShockTube.F90 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/GRHydro_ShockTube.F90 2011-04-27 22:25:23 UTC (rev 123) @@ -93,6 +93,17 @@ velzr = 0.d0 epsl = 1.5d0 epsr = 1.2d0 + else if (CCTK_EQUALS(shock_case,"Balsaralike1")) then + rhol = 1.0d0 + rhor = 0.125d0 + velxl = 0.0d0 + velxr = 0.0d0 + velyl = 0.0d0 + velyr = 0.0d0 + velzl = 0.0d0 + velzr = 0.0d0 + epsl = 1.0d0/rhol + epsr = 0.1d0/rhor !!$This line only for polytrope, k=1 !!$ epsr = 0.375d0 else if (CCTK_EQUALS(shock_case,"Blast")) then File [modified]: GRHydro_ShockTubeM.F90 Delta lines: +565 -8 =================================================================== --- trunk/src/GRHydro_ShockTubeM.F90 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/GRHydro_ShockTubeM.F90 2011-04-27 22:25:23 UTC (rev 123) @@ -23,6 +23,11 @@ #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) +#define OOSQRT2 (0.7071067811865475244008442) +#define OOSQRT3 (0.5773502691896257645091489) +#define OOSQRT6 (0.4082482904638630163662140) + + /*@@ @routine GRHydro_shocktubeM @date Sat Jan 26 02:53:49 2002 @@ -53,7 +58,9 @@ CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, & velzl, velzr, epsl, epsr CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr + CCTK_REAL :: ux,uy,uz,ut,tmp,tmp2,tmp3 + bvcxl = Bx_init bvcyl = By_init bvczl = Bz_init @@ -69,18 +76,27 @@ do j=1,ny do k=1,nz if (CCTK_EQUALS(shocktube_type,"diagshock")) then - direction = x(i,j,k) - shock_xpos + & - y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos +!!$ The diagshock choice yields a shock plane perpendicular to the fixed vector (1,1,1) +!!$ This could be changed, but would require 3 new params containing the new shock direction + direction = x(i,j,k) - shock_xpos + & + y(i,j,k) - shock_ypos + z(i,j,k) - shock_zpos + + else if (CCTK_EQUALS(shocktube_type,"diagshock2d")) then +!!$ The diagshock choice yields a shock plane perpendicular to the fixed vector (1,1,0), with similarity in the z-dir. +!!$ This could be changed, but would require 2 new params containing the new shock direction + direction = x(i,j,k) - shock_xpos + & + y(i,j,k) - shock_ypos + else if (CCTK_EQUALS(shocktube_type,"xshock")) then - direction = x(i,j,k) - shock_xpos + direction = x(i,j,k) - shock_xpos else if (CCTK_EQUALS(shocktube_type,"yshock")) then - direction = y(i,j,k) - shock_ypos + direction = y(i,j,k) - shock_ypos else if (CCTK_EQUALS(shocktube_type,"zshock")) then - direction = z(i,j,k) - shock_zpos + direction = z(i,j,k) - shock_zpos else if (CCTK_EQUALS(shocktube_type,"sphere")) then - direction = sqrt((x(i,j,k)-shock_xpos)**2+& - (y(i,j,k)-shock_ypos)**2+& - (z(i,j,k)-shock_zpos)**2)-shock_radius + direction = sqrt((x(i,j,k)-shock_xpos)**2+& + (y(i,j,k)-shock_ypos)**2+& + (z(i,j,k)-shock_zpos)**2)-shock_radius end if if (CCTK_EQUALS(shock_case,"Simple")) then rhol = 10.d0 @@ -117,6 +133,31 @@ velzr = 0.d0 epsl = 1500.d0 epsr = 1.5d-2 + +!!$ The following shocktubes are from Balsara 2001 . +!!$ All use n=1600 cells, over domain x=[-0.5,0.5] +!!$ All assume ideal-gas or gamma-law EOS, the first test uses GAMMA=2. while the rest GAMMA=5./3. + +!!$ Unmagnetized Test 1 (rel. Brio & Wu 1988 by van Putten 1993) of Balsara 2001 -- compare at t=0.4 + else if (CCTK_EQUALS(shock_case,"Balsara0")) then + rhol = 1.0d0 + rhor = 0.125d0 + velxl = 0.0d0 + velxr = 0.0d0 + velyl = 0.0d0 + velyr = 0.0d0 + velzl = 0.0d0 + velzr = 0.0d0 + bvcxl=0.0d0 + bvcxr=0.0d0 + bvcyl=0.0d0 + bvcyr=0.0d0 + bvczl=0.0d0 + bvczr=0.0d0 + epsl = 1.0d0/rhol + epsr = 0.1d0/rhor + +!!$ Test 1 (rel. Brio & Wu 1988 by van Putten 1993) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara1")) then rhol = 1.0d0 rhor = 0.125d0 @@ -134,6 +175,8 @@ bvczr=0.0d0 epsl = 1.0d0/rhol epsr = 0.1d0/rhor + +!!$ Test 2 (blast wave) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara2")) then rhol = 1.d0 rhor = 1.d0 @@ -151,6 +194,8 @@ bvczr=0.7d0 epsl = 1.5d0*30.0d0/rhol epsr = 1.5d0*1.0d0/rhor + +!!$ Test 3 (blast wave) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara3")) then rhol = 1.d0 rhor = 1.d0 @@ -168,6 +213,8 @@ bvczr=0.7d0 epsl = 1.5d0*1000.0d0/rhol epsr = 1.5d0*0.1d0/rhor + +!!$ Test 4 (rel. version of Noh 1987) of Balsara 2001 -- compare at t=0.4 else if (CCTK_EQUALS(shock_case,"Balsara4")) then rhol = 1.d0 rhor = 1.d0 @@ -185,6 +232,8 @@ bvczr=-7.d0 epsl = 1.5d0*0.1d0/rhol epsr = 1.5d0*0.1d0/rhor + +!!$ Test 5 (non-coplanar set of waves) of Balsara 2001 -- compare at t=0.55 else if (CCTK_EQUALS(shock_case,"Balsara5")) then rhol = 1.08d0 rhor = 1.d0 @@ -202,12 +251,286 @@ bvczr=0.5d0 epsl = 1.5d0*0.95d0/rhol epsr = 1.5d0*1.0d0/rhor + +!!$ "Generic Alfven Test of Giacomazzo and Rezzolla J.Comp.Phys (2006) + else if (CCTK_EQUALS(shock_case,"Alfven")) then + rhol = 1.d0 + rhor = 0.9d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.3d0 + velyr = 0.d0 + velzl = 0.4d0 + velzr = 0.d0 + bvcxl=1.0d0 + bvcxr=1.0d0 + bvcyl=6.d0 + bvcyr=5.d0 + bvczl=2.d0 + bvczr=2.d0 + epsl = 1.5d0*5.d0/rhol + epsr = 1.5d0*5.3d0/rhor + +!!$ The following 9 tests are from Komissarov 1999 . +!!$ Note that the data is specified in terms of the 4-velocity, so some conversion is necessary. +!!$ All assume ideal-gas or gamma-law EOS with GAMMA=4./3. + +!!$ Fast Shock test of Komissarov 1999 -- compare at t=2.5 , n=40, x=[-1,1] + else if (CCTK_EQUALS(shock_case,"Komissarov1")) then + rhol = 1.d0 + rhor = 25.48d0 + bvcxl=20.d0 + bvcxr=20.d0 + bvcyl=25.02d0 + bvcyr=49.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 367.5d0 /rhor + + ux=25.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=1.091d0 + uy=0.3923d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Slow Shock test of Komissarov 1999 -- compare at t=2. , n=200, x=[-0.5,1.5] + else if (CCTK_EQUALS(shock_case,"Komissarov2")) then + rhol = 1.d0 + rhor = 3.323d0 + bvcxl=10.d0 + bvcxr=10.d0 + bvcyl=18.28d0 + bvcyr=14.49d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 55.36d0 /rhor + + ux=1.53d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=0.9571d0 + uy=-0.6822d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Switch-off Fast test of Komissarov 1999 -- compare at t=1. , n=150, x=[-1,1] + else if (CCTK_EQUALS(shock_case,"Komissarov3")) then + rhol = 0.1d0 + rhor = 0.562d0 + bvcxl=2.d0 + bvcxr=2.d0 + bvcyl=0.d0 + bvcyr=4.710d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 10.d0 /rhor + + ux=-2.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=-0.212d0 + uy=-0.590d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Switch-on Slow test of Komissarov 1999 -- compare at t=2. , n=150, x=[-1,1.5] + else if (CCTK_EQUALS(shock_case,"Komissarov4")) then + rhol = 1.78d-2 + rhor = 1.d-2 + bvcxl=1.d0 + bvcxr=1.d0 + bvcyl=1.022d0 + bvcyr=0.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 0.1d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=-0.765d0 + uy=-1.386d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=0.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Alfven wave test of Komissarov 1999 -- compare at t=2. , n=200, x=[-1,1.5] +!!$ Needs special setup -- FIX + else if (CCTK_EQUALS(shock_case,"Komissarov5")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=3.d0 + bvcxr=3.d0 + bvcyl=3.d0 + bvcyr=-6.857d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=0.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=3.70d0 + uy=5.76d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Compound wave test of Komissarov 1999 -- compare at t=0.1,0.75,1.5 , n=200, x=[-0.5,1.5] +!!$ Needs special setup -- FIX + else if (CCTK_EQUALS(shock_case,"Komissarov6")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=3.d0 + bvcxr=3.d0 + bvcyl=3.d0 + bvcyr=-6.857d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=0.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=3.70d0 + uy=5.76d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + +!!$ Shock Tube 1 test of Komissarov 1999 -- compare at t=1. , n=400, x=[-1,1.5] + else if (CCTK_EQUALS(shock_case,"Komissarov7")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=1.d0 + bvcxr=1.d0 + bvcyl=0.d0 + bvcyr=0.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d3 /rhol + epsr = 3.d0 * 1.d0 /rhor + + velxl = 0. + velyl = 0. + velzl = 0. + + velxr = 0. + velyr = 0. + velzr = 0. + +!!$ Shock Tube 2 test of Komissarov 1999 -- compare at t=1. , n=500, x=[-1.25,1.25] + else if (CCTK_EQUALS(shock_case,"Komissarov8")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=0.d0 + bvcxr=0.d0 + bvcyl=20.d0 + bvcyr=0.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 30.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + velxl = 0. + velyl = 0. + velzl = 0. + + velxr = 0. + velyr = 0. + velzr = 0. + +!!$ Collision test of Komissarov 1999 -- compare at t=1.22 , n=200, x=[-1,1] + else if (CCTK_EQUALS(shock_case,"Komissarov9")) then + rhol = 1.d0 + rhor = 1.d0 + bvcxl=10.d0 + bvcxr=10.d0 + bvcyl=10.d0 + bvcyr=-10.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 3.d0 * 1.d0 /rhol + epsr = 3.d0 * 1.d0 /rhor + + ux=5.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxl = ux/ut + velyl = uy/ut + velzl = uz/ut + + ux=-5.d0 + uy=0.d0 + uz=0.d0 + ut = sqrt(1. + ux*ux + uy*uy + uz*uz) + velxr = ux/ut + velyr = uy/ut + velzr = uz/ut + else call CCTK_WARN(0,"Shock case not recognized") end if if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then + +!!$ Left state + rho(i,j,k) = rhol velx(i,j,k) = velxl vely(i,j,k) = velyl @@ -217,6 +540,9 @@ Bvecy(i,j,k)=bvcyl Bvecz(i,j,k)=bvczl else + +!!$ Right state + rho(i,j,k) = rhor velx(i,j,k) = velxr vely(i,j,k) = velyr @@ -227,6 +553,58 @@ Bvecz(i,j,k)=bvczr end if + + if (CCTK_EQUALS(shocktube_type,"yshock")) then +!!$ Cycle x,y,z forward + tmp=velx(i,j,k) + velx(i,j,k)=velz(i,j,k) + velz(i,j,k)=vely(i,j,k) + vely(i,j,k)=tmp + tmp=Bvecx(i,j,k) + Bvecx(i,j,k)=Bvecz(i,j,k) + Bvecz(i,j,k)=Bvecy(i,j,k) + Bvecy(i,j,k)=tmp + else if (CCTK_EQUALS(shocktube_type,"zshock")) then +!!$ Cycle x,y,z backward + tmp=velx(i,j,k) + velx(i,j,k)=vely(i,j,k) + vely(i,j,k)=velz(i,j,k) + velz(i,j,k)=tmp + tmp=Bvecx(i,j,k) + Bvecx(i,j,k)=Bvecy(i,j,k) + Bvecy(i,j,k)=Bvecz(i,j,k) + Bvecz(i,j,k)=tmp + else if (CCTK_EQUALS(shocktube_type,"diagshock")) then +!!$ Rotated basis vectors necessary to evaluate the orthogonal matrix elements: +!!$ xhat = 1/sqrt(3)[1,1,1], yhat = 1/sqrt(2)[-1,1,0]; zhat = 1/sqrt(6)[-1,-1,2] +!!$ Orthogonal matrix constructed from the tensor product between the original +!!$ cartesian basis x's and new basis vectors xhat's, rotated towards the diagonal +!!$ shock normal and tangent directions. + tmp = OOSQRT3*velx(i,j,k) - OOSQRT2*vely(i,j,k) - OOSQRT6*velz(i,j,k) + tmp2 = OOSQRT3*velx(i,j,k) + OOSQRT2*vely(i,j,k) - OOSQRT6*velz(i,j,k) + tmp3 = OOSQRT3*velx(i,j,k) + 2.d0*OOSQRT6*velz(i,j,k) + velx(i,j,k)=tmp + vely(i,j,k)=tmp2 + velz(i,j,k)=tmp3 + tmp = OOSQRT3*Bvecx(i,j,k) - OOSQRT2*Bvecy(i,j,k) - OOSQRT6*Bvecz(i,j,k) + tmp2 = OOSQRT3*Bvecx(i,j,k) + OOSQRT2*Bvecy(i,j,k) - OOSQRT6*Bvecz(i,j,k) + tmp3 = OOSQRT3*Bvecx(i,j,k) + 2.d0*OOSQRT6*Bvecz(i,j,k) + Bvecx(i,j,k)=tmp + Bvecy(i,j,k)=tmp2 + Bvecz(i,j,k)=tmp3 + else if (CCTK_EQUALS(shocktube_type,"diagshock2d")) then +!!$ New basis: +!!$ xhat = 1/sqrt(2)[1,1,0], yhat = 1/sqrt(2)[-1,1,0]; zhat = [0,0,1] + tmp = OOSQRT2*velx(i,j,k) - OOSQRT2*vely(i,j,k) + tmp2 = OOSQRT2*velx(i,j,k) + OOSQRT2*vely(i,j,k) + velx(i,j,k)=tmp + vely(i,j,k)=tmp2 + tmp = OOSQRT2*Bvecx(i,j,k) - OOSQRT2*Bvecy(i,j,k) + tmp2 = OOSQRT2*Bvecx(i,j,k) + OOSQRT2*Bvecy(i,j,k) + Bvecx(i,j,k)=tmp + Bvecy(i,j,k)=tmp2 + endif + 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)) if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then @@ -257,3 +635,182 @@ return end subroutine GRHydro_shocktubeM + +subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew + CCTK_INT :: xoff,yoff,zoff,indsum + CCTK_REAL :: det + + stenp1=GRHydro_stencil + 1 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + minsum = 3*stenp1 + maxsum = nx+ny+nz-3*stenp1 + 3 + + xoff=0 + yoff=0 + zoff=0 + + do k=1,nz + if(k.lt.stenp1)then + zoff=k-stenp1 + else if(k.gt.nz-stenp1+1) then + zoff=k-(nz-stenp1+1) + else + zoff=0 + endif + + do j=1,ny + if(j.lt.stenp1)then + yoff=j-stenp1 + else if(j.gt.ny-stenp1+1) then + yoff=j-(ny-stenp1+1) + else + yoff=0 + endif + + do i=1,nx + if(i.lt.stenp1) then + xoff=i-stenp1 + else if(i.gt.nx-stenp1+1) then + xoff=i-(nx-stenp1+1) + else + xoff=0 + endif + + indsum = i+j+k + + if( (xoff.ne.0.or.yoff.ne.0.or.zoff.ne.0) .and. & + indsum.ge.minsum.and.indsum.le.maxsum) then +!!$ We can map the point to the interior diagonal, orthogonal to the shock. + + inew=indsum/3 + jnew=(indsum-inew)/2 + knew=indsum-inew-jnew + + dens(i,j,k) = dens(inew,jnew,knew) + sx(i,j,k) = sx(inew,jnew,knew) + sy(i,j,k) = sy(inew,jnew,knew) + sz(i,j,k) = sz(inew,jnew,knew) + tau(i,j,k) = tau(inew,jnew,knew) + Bvecx(i,j,k)=Bvecx(inew,jnew,knew) + Bvecy(i,j,k)=Bvecy(inew,jnew,knew) + Bvecz(i,j,k)=Bvecz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + endif + + enddo + enddo + enddo + +end subroutine GRHydro_Diagshock_BoundaryM + +subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, kc, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew + CCTK_INT :: xoff,yoff,zoff,indsum + CCTK_REAL :: det + + stenp1=GRHydro_stencil + 1 + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + minsum = 2*stenp1 + maxsum = nx+ny-2*stenp1 + 2 + + xoff=0 + yoff=0 + zoff=0 + + do k=1,nz + if(k.lt.stenp1)then + zoff=k-stenp1 + else if(k.gt.nz-stenp1+1) then + zoff=k-(nz-stenp1+1) + else + zoff=0 + endif + + do j=1,ny + if(j.lt.stenp1)then + yoff=j-stenp1 + else if(j.gt.ny-stenp1+1) then + yoff=j-(ny-stenp1+1) + else + yoff=0 + endif + + do i=1,nx + if(i.lt.stenp1) then + xoff=i-stenp1 + else if(i.gt.nx-stenp1+1) then + xoff=i-(nx-stenp1+1) + else + xoff=0 + endif + + indsum = i+j + + if( (xoff.ne.0.or.yoff.ne.0) .and. & + indsum.ge.minsum.and.indsum.le.maxsum) then +!!$ We can map the point to the interior diagonal, orthogonal to the shock. + + inew=indsum/2 + jnew=indsum-inew + + dens(i,j,k) = dens(inew,jnew,k) + sx(i,j,k) = sx(inew,jnew,k) + sy(i,j,k) = sy(inew,jnew,k) + sz(i,j,k) = sz(inew,jnew,k) + tau(i,j,k) = tau(inew,jnew,k) + Bvecx(i,j,k)=Bvecx(inew,jnew,k) + Bvecy(i,j,k)=Bvecy(inew,jnew,k) + Bvecz(i,j,k)=Bvecz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + else if( zoff.ne.0) then + + kc = nz/2 + + dens(i,j,k) = dens(i,j,kc) + sx(i,j,k) = sx(i,j,kc) + sy(i,j,k) = sy(i,j,kc) + sz(i,j,k) = sz(i,j,kc) + tau(i,j,k) = tau(i,j,kc) + Bvecx(i,j,k)=Bvecx(i,j,kc) + Bvecy(i,j,k)=Bvecy(i,j,kc) + Bvecz(i,j,k)=Bvecz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + endif + + enddo + enddo + enddo + +end subroutine GRHydro_Diagshock2D_BoundaryM File [modified]: make.code.defn Delta lines: +3 -1 =================================================================== --- trunk/src/make.code.defn 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/src/make.code.defn 2011-04-27 22:25:23 UTC (rev 123) @@ -14,7 +14,9 @@ CheckParam.c\ GRHydro_P2C2P.F90 \ GRHydro_P2C2PM.F90 \ - GRHydro_P2C2PM_polytype.F90 + GRHydro_P2C2PM_polytype.F90 \ + GRHydro_MonopoleM.F90 \ + GRHydro_CylindricalExplosionM.F90 # Subdirectories containing source files SUBDIRS = Directory: /trunk/ ================== File [modified]: param.ccl Delta lines: +18 -1 =================================================================== --- trunk/param.ccl 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/param.ccl 2011-04-27 22:25:23 UTC (rev 123) @@ -12,6 +12,8 @@ "only_atmo" :: "Set only a low atmosphere" "read_conformal":: "After reading in initial alp, rho and gxx from h5 files, sets the other quantities" "simple_wave" :: "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" + "Monopole" :: "Monopole at the center" + "cylexp" :: "Cylindrical Explosion" } shares:ADMBase @@ -31,22 +33,35 @@ KEYWORD shocktube_type "Diagonal or parallel shock?" { "diagshock" :: "Diagonal across all axes" + "diagshock2d" :: "Diagonal across x-y axes" "xshock" :: "Parallel to x axis" "yshock" :: "Parallel to y axis" "zshock" :: "Parallel to z axis" "sphere" :: "spherically symmetric shock" -} "diagshock" +} "xshock" KEYWORD shock_case "Simple, Sod's problem or other?" { "Simple" :: "GRAstro_Hydro test case" "Sod" :: "Sod's problem" "Blast" :: "Strong blast wave" + "Balsaralike1" :: "Hydro version of Balsara Test #1" + "Balsara0" :: "Balsara Test #1, but unmagnetized" "Balsara1" :: "Balsara Test #1" "Balsara2" :: "Balsara Test #2" "Balsara3" :: "Balsara Test #3" "Balsara4" :: "Balsara Test #4" "Balsara5" :: "Balsara Test #5" + "Alfven" :: "Generical Alfven Test" + "Komissarov1" :: "Komissarov Test #1" + "Komissarov2" :: "Komissarov Test #2" + "Komissarov3" :: "Komissarov Test #3" + "Komissarov4" :: "Komissarov Test #4" + "Komissarov5" :: "Komissarov Test #5" + "Komissarov6" :: "Komissarov Test #6" + "Komissarov7" :: "Komissarov Test #7" + "Komissarov8" :: "Komissarov Test #8" + "Komissarov9" :: "Komissarov Test #9" } "Sod" REAL shock_xpos "Position of shock plane: x" @@ -124,3 +139,5 @@ USES REAL initial_rho_abs_min USES REAL initial_rho_rel_min USES REAL initial_atmosphere_factor +USES int GRHydro_stencil +USES BOOLEAN clean_divergence File [modified]: schedule.ccl Delta lines: +30 -0 =================================================================== --- trunk/schedule.ccl 2011-04-21 04:46:47 UTC (rev 122) +++ trunk/schedule.ccl 2011-04-27 22:25:23 UTC (rev 123) @@ -5,6 +5,13 @@ LANG: C } "Check parameters" +if (CCTK_Equals(initial_hydro,"Monopole")) { + schedule GRHydro_MonopoleM in HydroBase_Initial + { + LANG: Fortran + } "Monopole initial data - MHD version" +} + if (CCTK_Equals(initial_hydro,"shocktube")) { if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) { @@ -13,6 +20,22 @@ LANG: Fortran } "Shocktube initial data - MHD version" + if(CCTK_Equals(shocktube_type,"diagshock")) + { + schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + } "Diagonal shock boundary conditions" + } + + if(CCTK_Equals(shocktube_type,"diagshock2d")) + { + schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + } "2-D Diagonal shock boundary conditions" + } + } else { schedule GRHydro_shocktube in HydroBase_Initial { @@ -22,6 +45,13 @@ } } +if (CCTK_Equals(initial_hydro,"cylexp")) { + schedule GRHydro_CylindricalExplosionM in HydroBase_Initial + { + LANG: Fortran + } "Cylindrical Explosion initial data - MHD-only" +} + if (CCTK_Equals(initial_data,"con2primtest")) { STORAGE:GRHydro_init_data_reflevel schedule GRHydro_Init_Data_RefinementLevel IN HydroBase_Initial BEFORE GRHydro_con2primtest

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 233)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 05:20 PM Modified: /branches/divclean/src/ GRHydro_Con2Prim.F90, GRHydro_Con2PrimM.F90 Log: Fix wrongly overwritten Con2Prim. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_Con2Prim.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2PrimM.F90 Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/ (Rev. 232)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 04:02 PM Modified: /branches/divclean/ interface.ccl, schedule.ccl /branches/divclean/src/ GRHydro_Boundaries.F90, GRHydro_CalcUpdate.F90, GRHydro_Con2Prim.F90, GRHydro_HLLEM.F90, GRHydro_Prim2ConM.F90, GRHydro_Reconstruct.F90, GRHydro_ReconstructPoly.F90, GRHydro_RegisterVars.cc, GRHydro_UpdateMaskM.F90 Log: Add Josh's changes. File Changes: Directory: /branches/divclean/ ============================== File [modified]: interface.ccl Delta lines: <Unavailable> File [modified]: schedule.ccl Delta lines: <Unavailable> Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_Boundaries.F90 Delta lines: <Unavailable> File [modified]: GRHydro_CalcUpdate.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Con2Prim.F90 Delta lines: <Unavailable> File [modified]: GRHydro_HLLEM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Prim2ConM.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Reconstruct.F90 Delta lines: <Unavailable> File [modified]: GRHydro_ReconstructPoly.F90 Delta lines: <Unavailable> File [modified]: GRHydro_RegisterVars.cc Delta lines: <Unavailable> File [modified]: GRHydro_UpdateMaskM.F90 Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/branches/divclean/src/ (Rev. 231)
by tanja.bode＠physics.gatech.edu 27 Apr '11

27 Apr '11

User: tbode Date: 2011/04/27 03:30 PM Modified: /branches/divclean/src/ GRHydro_SourceM.F90 Log: Modify rhs for divergence cleaning. File Changes: Directory: /branches/divclean/src/ ================================== File [modified]: GRHydro_SourceM.F90 Delta lines: +50 -3 =================================================================== --- branches/divclean/src/GRHydro_SourceM.F90 2011-04-27 17:51:47 UTC (rev 230) +++ branches/divclean/src/GRHydro_SourceM.F90 2011-04-27 20:30:04 UTC (rev 231) @@ -33,6 +33,12 @@ #define Bvecx(i,j,k) Bvec(i,j,k,1) #define Bvecy(i,j,k) Bvec(i,j,k,2) #define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) +#define Bconsrhsx(i,j,k) Bconsrhs(i,j,k,1) +#define Bconsrhsy(i,j,k) Bconsrhs(i,j,k,2) +#define Bconsrhsz(i,j,k) Bconsrhs(i,j,k,3) /*@@ @routine SourceTermsM @@ -76,8 +82,9 @@ CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz CCTK_REAL :: invalp, invalp2 - CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc + !!CCTK_REAL :: dx_psidc, dy_psidc, dz_psidc CCTK_REAL :: bvcx_source, bvcy_source, bvcz_source + CCTK_REAL :: dx_det_bydet, dy_det_bydet, dz_det_bydet CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow,bdotv,b2,dum,bxlow,bylow,bzlow CCTK_REAL :: bt,bx,by,bz,rhohstarW2,pstar @@ -103,7 +110,7 @@ densrhs = 0.d0 srhs = 0.d0 taurhs = 0.d0 - Bvecrhs=0.d0 + Bconsrhs=0.d0 if (evolve_tracer .ne. 0) then cons_tracerrhs = 0.d0 @@ -436,7 +443,47 @@ srhs(i,j,k,3) = alp(i,j,k)*sqrtdet*sz_source taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source if(clean_divergence.ne.0) then - psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k) + psidcrhs(i,j,k) = -1.d0 * ( ch_dc*ch_dc/cp_dc/cp_dc*alp(i,j,k) + & + dx_shiftx + dy_shifty + dz_shiftz ) * psidc(i,j,k) + & + Bconsx(i,j,k) * ( dx_alp - half*alp(i,j,k) * & + ( uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + 2.d0*uxy*dx_gxy + & + 2.d0*uxz*dx_gxz + 2.d0*uyz*dx_gyz ) )/ sqrtdet + & + Bconsy(i,j,k) * (dy_alp - half*alp(i,j,k) * & + ( uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & + 2.d0*uxz*dy_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + & + Bconsz(i,j,k) * (dz_alp - half*alp(i,j,k) * & + ( uxx*dz_gxx + uyy*dy_gyy + uzz*dy_gzz + 2.d0*uxy*dy_gxy + & + 2.d0*uxz*dz_gxz + 2.d0*uyz*dy_gyz ) )/ sqrtdet + + !!$ g^{jk} d_i g_{kj} = d_i (g) / det + dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + & + 2.d0*(uxy*dx_gxy+uxz*dx_gxz+uyz*dx_gyz) + dy_det_bydet = uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + & + 2.d0*(uxy*dy_gxy+uxz*dy_gxz+uyz*dy_gyz) + dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + & + 2.d0*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz) + + bvcx_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftx + & + Bconsy(i,j,k)*dy_shiftx + Bconsz(i,j,k)*dz_shiftx ) + & + psidc(i,j,k)*sqrtdet*( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + & + psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxx*dx_det_bydet + & + uxy*dy_det_bydet + uxz*dz_det_bydet ) + + bvcy_source = -1.d0 * ( Bconsx(i,j,k)*dx_shifty + & + Bconsy(i,j,k)*dy_shifty + Bconsz(i,j,k)*dz_shifty ) + & + psidc(i,j,k)*sqrtdet*( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + & + psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxy*dx_det_bydet + & + uyy*dy_det_bydet + uyz*dz_det_bydet ) + + bvcz_source = -1.d0 * ( Bconsx(i,j,k)*dx_shiftz + & + Bconsy(i,j,k)*dy_shiftz + Bconsz(i,j,k)*dz_shiftz ) + & + psidc(i,j,k)*sqrtdet*( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + & + psidc(i,j,k)*3.d0*half*alp(i,j,k)*sqrtdet*( uxz*dx_det_bydet + & + uyz*dy_det_bydet + uzz*dz_det_bydet ) + + Bconsrhsx(i,j,k) = bvcx_source + Bconsrhsy(i,j,k) = bvcy_source + Bconsrhsz(i,j,k) = bvcz_source endif enddo

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 71)
by cott＠tapir.caltech.edu 27 Apr '11

27 Apr '11

User: cott Date: 2011/04/27 01:57 PM Modified: / ET.tex /local_bibtex/ ott_references.bib Log: * add a ref File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +1 -1 =================================================================== --- ET.tex 2011-04-27 18:54:54 UTC (rev 70) +++ ET.tex 2011-04-27 18:57:16 UTC (rev 71) @@ -965,7 +965,7 @@ specific internal energy and the pressure, respectively. \subsubsection{Spacetime Curvature Evolution} The Einstein Toolkit -curvature evolution code \codename{McLachlan}~\cite{brown:09} is +curvature evolution code \codename{McLachlan}~\cite{brown:09,reisswig:11ccwave} is auto-generated from tensor equations via \codename{Kranc} (Sec.~\ref{sec:kranc}). It implements the Einstein equations in a $3+1$ split as a Cauchy initial boundary value problem Directory: /local_bibtex/ ========================= File [modified]: ott_references.bib Delta lines: +12 -0 =================================================================== --- local_bibtex/ott_references.bib 2011-04-27 18:54:54 UTC (rev 70) +++ local_bibtex/ott_references.bib 2011-04-27 18:57:16 UTC (rev 71) @@ -165,3 +165,15 @@ fulltexturl = {http://link.aps.org/abstract/PRD/v79/e044023}, doi = {10.1103/PhysRevD.79.044023}, } + +@ARTICLE{reisswig:11ccwave, + author = {{Reisswig}, C. and {Ott}, C.~D. and {Sperhake}, U. and {Schnetter}, E.}, + title = "{Gravitational wave extraction in simulations of rotating stellar core collapse}", + journal = {\prd}, + keywords = {Numerical relativity, Finite-difference methods, Wave generation and sources, Supernovae}, + year = 2011, + month = mar, + volume = 83, + number = 6, + pages = {064008}, +}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 70)
by cott＠tapir.caltech.edu 27 Apr '11

27 Apr '11

User: cott Date: 2011/04/27 01:54 PM Modified: / ET.tex /local_bibtex/ ott_references.bib Log: * more text on McLachlan File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +27 -25 =================================================================== --- ET.tex 2011-04-27 18:44:41 UTC (rev 69) +++ ET.tex 2011-04-27 18:54:54 UTC (rev 70) @@ -1157,18 +1157,18 @@ \subsubsection{Boundary Conditions} -During time evolution, we apply a Sommerfeld-type radiative boundary -condition to all components of the evolved BSSN variables as described -in \cite{Alcubierre2000}. The main feature of this boundary condition -is that it assumes approximate spherical symmetry of the solution, -while applying the actual boundary condition on the boundary of a -cubic grid where the face normals are not aligned with the radial -direction. This boundary condition defines the right hand side -of the BSSN state vector on the outer boundary, which is then -integrated in time as well, so that the boundary and interior are -calculated with the same order of accuracy. +During time evolution, a Sommerfeld-type radiative boundary condition +is applied to all components of the evolved BSSN variables as +described in \cite{Alcubierre2000}. The main feature of this boundary +condition is that it assumes approximate spherical symmetry of the +solution, while applying the actual boundary condition on the boundary +of a cubic grid where the face normals are not aligned with the radial +direction. This boundary condition defines the right hand side of the +BSSN state vector on the outer boundary, which is then integrated in +time as well, so that the boundary and interior are calculated with +the same order of accuracy. -The main part of the boundary condition assumes that we have an +The main part of the boundary condition assumes that one has an outgoing radial wave with some speed $v_0$: \begin{eqnarray} X & = & X_0 + \frac{u(r - v_0 t)}{r} @@ -1182,17 +1182,18 @@ \end{eqnarray} where $v^i = v_0\, x^i/r$. The spatial derivatives $\partial_i$ are evaluated using centered finite differencing where possible, and -one-sided finite differencing elsewhere. We use second order stencils -in our implementation. +one-sided finite differencing elsewhere. Second order stencils +are used in the current implementation. -In addition to this main part, we also account for those parts of the -solution that do not behave as a pure wave, e.g., Coulomb type terms -caused by infall of the coordinate lines. We assume that these parts -decay with a certain power $p$ of the radius. We implement this by -considering the radial derivative of the source term above, and -extrapolating according to this power-law decay. +In addition to this main part, it is also necessary to account for +those parts of the solution that do not behave as a pure wave, e.g., +Coulomb type terms caused by infall of the coordinate lines. The +assumption is made that these parts decay with a certain power $p$ of +the radius. This is implemented by considering the radial derivative of +the source term above, and extrapolating according to this power-law +decay. -Given a source term $(\partial_t X)$, we define the corrected source +Given a source term $(\partial_t X)$, one defines the corrected source term $(\partial_t X)^*$ via \begin{eqnarray} (\partial_t X)^* & = & (\partial_t X) + \left( \frac{r}{r - n^i @@ -1201,8 +1202,8 @@ where $n^i$ is the normal vector of the corresponding boundary face. The spatial derivatives $\partial_i$ are evaluated by comparing neighbouring grid points, corresponding to a second-order stencil -evaluated in the middle between the two neighbouring grid points. We -assume a second-order decay, i.e., we choose $p=2$. +evaluated in the middle between the two neighbouring grid points. +Second-order decays is assumed, hence $p=2$. As with the initial conditions above, this boundary condition is evaluated on several layers of grid points, starting from the @@ -1212,19 +1213,20 @@ This boundary condition is only a coarse approximation of the actual decay behavior of the BSSN state vector, and it does not capture the -correct behavior of the evolved variables. However, we observe that +correct behavior of the evolved variables. However, one finds that this boundary condition leads to stable evolutions if applied sufficiently far from the source. Errors introduced at the boundary (both errors in the geometry and constraint violations) propagate inwards with the speed of light \cite{ES-Brown2007b}. Gauge changes introduced by the boundary condition, which are physically not -observable, propagate faster, with a speed up to $\sqrt{2}$ for our -gauge conditions. +observable, propagate faster, with a speed up to $\sqrt{2}$ for the + gauge conditions used in \codename{McLachlan}. + \subsubsection{Hydrodynamics: \codename{GRHydro}} \label{sec:GRHydro} Directory: /local_bibtex/ ========================= File [modified]: ott_references.bib Delta lines: +25 -0 =================================================================== --- local_bibtex/ott_references.bib 2011-04-27 18:44:41 UTC (rev 69) +++ local_bibtex/ott_references.bib 2011-04-27 18:54:54 UTC (rev 70) @@ -140,3 +140,28 @@ doi = {10.1088/0264-9381/27/16/167001}, } + +@article{Alcubierre2000, + author = "Alcubierre, M. and Br{\"u}gmann, B. and Dramlitsch, T. and Font, J. A. and Papadopoulos, P. and Seidel, E. and Stergioulas, N. and Takahashi, R.", + title = "Towards a Stable Numerical Evolution of Strongly Gravitating Systems in General Relativity: The Conformal Treatments", + journal = "Phys. Rev. D", + year = "2000", + volume = "62", + pages = "044034"} + + + +@Article{ES-Brown2007b, + status = {refereed}, + author = {David Brown and Peter Diener and Olivier Sarbach and + Erik Schnetter and Manuel Tiglio}, + title = {Turduckening black holes: an analytical and + computational study}, + journal = {Phys. Rev. D}, + year = 2009, + volume = 79, + pages = 044023, + receiveddate = {2008-09-20}, + fulltexturl = {http://link.aps.org/abstract/PRD/v79/e044023}, + doi = {10.1103/PhysRevD.79.044023}, +}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 69)
by diener＠cct.lsu.edu 27 Apr '11

27 Apr '11

User: diener Date: 2011/04/27 01:44 PM Modified: / ET.tex Log: Some stuff in relativity tools. File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +25 -5 =================================================================== --- ET.tex 2011-04-27 18:39:57 UTC (rev 68) +++ ET.tex 2011-04-27 18:44:41 UTC (rev 69) @@ -1543,11 +1543,31 @@ \subsection{Relativity Tools\pages{3 Peter}} -\paragraph{Black Hole Excision} +\subsubsection*{Vacuum-Matter Coupling} +The thorn \codename{TmunuBase} provides grid functions for the stress-energy +tensor $T_{\mu\nu}$ as well as schedule groups to manage when $T_{\mu\nu}$ is +calculated. In a simulation many different thorns may contribute to the +stress-energy tensor and this thorn allows them to do so without explicit +dependence on each other. The resulting stress-energy tensor can then be used +by the spacetime evolution thorn (again without explicit dependence on any +matter thorns). When thorn \codename{MoL} is used for time evolution this +provides a high order spacetime-matter coupling. -\paragraph{Vacuum-Matter Coupling} +The grid functions provided by \codename{TmunuBase} are: +\begin{itemize} + \item The time compenent $T_{00}$ {\tt eTtt} + \item The mixed components $T_{0i}$ {\tt eTtx}, {\tt eTty}, {\tt eTtz} + \item The spatial components $T_{ij}$ {\tt eTxx}, {\tt eTxy}, {\tt eTxz}, {\tt eTyy}, {\tt eTyz}, {\tt eTzz} +\end{itemize} +In addition the grid scalar {\tt stress\_energy\_state} has the value 1 if storage is on for the stress-energy tensor and 0 if not. -\paragraph{Object tracking} + +Thorn \codename{ADMCoupling} provides a similar (but older) interface between +space-time and matter thorns. However, since it is based on an include file +mechansim it is more complicated to use and we recommend all new matter thorns +to use \codename{TmunuBase} instead. + +\subsubsection*{Object tracking} \label{sec:object-tracking} We provide a module (\codename{PunctureTracker}) for tracking black hole positions when evolved with the moving puncture technique. This can be used @@ -1565,8 +1585,8 @@ puncture location usually does not coincide with gridpoints, the shift is interpolated to the location of the puncture. Equation~(\ref{eq:puncturetracking}) is implemented with a simple first order -Euler scheme, which seems to be accurate enough for controlling the mesh -refinement grids. +Euler scheme, which seems to be accurate enough for controlling the location +of the mesh refinement hierarchy. \todo{ES: I think this paragraph provides too much detail for this paper.} \todo{ES: Can also track objects via their maximum density or centre of mass; I believe GRHydro does something like this}

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 68)
by cott＠tapir.caltech.edu 27 Apr '11

27 Apr '11

User: cott Date: 2011/04/27 01:39 PM Modified: / ET.tex Log: * more great stuff in curvature evolution section File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +11 -11 =================================================================== --- ET.tex 2011-04-27 18:24:03 UTC (rev 67) +++ ET.tex 2011-04-27 18:39:57 UTC (rev 68) @@ -696,6 +696,8 @@ a central access point for analysis thorns. \subsection{Initial Data\pages{4 Josh/Bruno}} +\label{sec:initial_data} + The Einstein Toolkit contains many modules used to generate initial data for general relativistic simulations, including both vacuum and hydrodynamical configurations. @@ -1131,16 +1133,14 @@ \subsubsection{Initial Conditions} -We set up our initial condition from the ADM variables $g_{ij}$, -$K_{ij}$, lapse $\alpha$, and shift $\beta^i$, as provided by the -initial data discussed in Sec.~\ref{sec:initial_models}. From these we -calculate the BSSN quantities via their definition, setting $B^i=0$, -and using cubic extrapolation for $\tilde\Gamma^i$ at the outer -boundary. This extrapolation is necessary since the $\tilde\Gamma^i$ are -calculated from derivatives of the metric, and one cannot use centered -finite differencing stencils near the outer boundary. We assume that -one could instead also use one-sided derivatives to calculate -$\tilde\Gamma^i$ on the boundary. +Initial conditions from the ADM variables $g_{ij}$, $K_{ij}$, lapse +$\alpha$, and shift $\beta^i$, as provided by the initial data +discussed in Sec.~\ref{sec:initial_data}. From these the BSSN +quantities are calculated via their definition, setting $B^i=0$, and +using cubic extrapolation for $\tilde\Gamma^i$ at the outer +boundary. This extrapolation is necessary since the $\tilde\Gamma^i$ +are calculated from derivatives of the metric, and one cannot use +centered finite differencing stencils near the outer boundary. The extrapolation stencils distinguish between points on the faces, edges, and corners of the grid. Points on the faces are extrapolated @@ -1150,7 +1150,7 @@ edge are extrapolated in the $(1,1,0)$ direction, while points in the $(+x,+y+z)$ corner are extrapolated in the $(1,1,1)$ direction. Since several layers of boundary points have to be filled for higher order -schemes (e.g., three layers for a fourth order scheme), we proceed +schemes (e.g., three layers for a fourth order scheme), one proceeds outwards starting from the innermost layer. Each subsequent layer is then defined via the points in the interior and the previously calculated layers.

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 67)
by cott＠tapir.caltech.edu 27 Apr '11

27 Apr '11

User: cott Date: 2011/04/27 01:24 PM Modified: / ET.tex /local_bibtex/ ott_references.bib Log: * clean up BSSN section a bit File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +12 -23 =================================================================== --- ET.tex 2011-04-27 17:58:18 UTC (rev 66) +++ ET.tex 2011-04-27 18:24:03 UTC (rev 67) @@ -931,11 +931,12 @@ In the following, we assume that the reader is familiar with the basics of numerical relativity and GR hydrodynamics. Detailed introductions to numerical relativity have recently been given by -Alcubierre~\cite{alcubierre:08}, Baumgarte \& Shapiro~\cite{baumgarte:10}, -and Centrella et al.~\cite{centrella:10}. +Alcubierre~\cite{alcubierre:08}, Baumgarte \& +Shapiro~\cite{baumgarte:10}, and Centrella et al.~\cite{centrella:10}. GR hydrodynamics has been reviewed by Font~\cite{font:08}. In the -following, we assume the reader to be familiar with general relativity. -We assume $G = c = M_\odot = 1$ throughout. +following, we assume the reader to be familiar with general +relativity, differential geometry and tensor analysis. We set $G = c += M_\odot = 1$ throughout. The Einstein Toolkit provides code to evolve the Einstein equations \begin{equation} @@ -957,8 +958,8 @@ T^{\mu\nu} = \rho h u^\mu u^\nu - g^{\mu\nu} P\,\,, \end{equation} where $\rho$ is the rest-mass density, $u^\mu$ is the 4-velocity, -$g^{\mu\nu}$ is the 4-metric, and $h = 1 + \epsilon + P/\rho$ is the -relativistic specific enthalpy with $\epsilon$ and $P$ being the +$g^{\mu\nu}$ is the upper 4-metric, and $h = 1 + \epsilon + P/\rho$ is +the relativistic specific enthalpy with $\epsilon$ and $P$ being the specific internal energy and the pressure, respectively. \subsubsection{Spacetime Curvature Evolution} The Einstein Toolkit @@ -971,7 +972,11 @@ \cite{shibata:95,baumgarte:95,alcubierre:00} of the original Arnowitt-Deser-Misner (ADM) formalism~\cite{adm:62} is employed. +\codename{McLachlan} uses fourth-order accurate finite differencing +for the spacetime variables and add a fifth-order Kreiss-Oliger +dissipation term to remove high frequency noise. + The evolved variables are the conformal factor $\Phi$, the conformal 3-metric $\tilde{\gamma}_{ij}$, the trace $K$ of the extrinsic curvature, the trace free extrinsic curvature $A_{ij}$ and the conformal connection @@ -1120,25 +1125,9 @@ with a constant $R$ defining the transition radius between the interior, where $q\approx1$, and the exterior, where $q$ falls off as $1/r$. Eq.~\ref{eq:eta} describes how $q$ appears in the gauge -parameters. In this paper we use $R=250\,M_\odot$ ($R = -369.2\,\mathrm{km}$). +parameters. -We implement the above BSSN equations and gauge conditions in the -{\tt McLachlan} code \cite{ES-Brown2007b, ES-mclachlanweb} which is -freely available as part of the {\tt EinsteinToolkit}. {\tt McLachlan} is -auto-generated from the definition of the variables and equations in the -{\tt Mathematica} format by the {\tt Kranc} code generator \cite{kranc04, - Husa:2004ip, krancweb}. {\tt Kranc} is a suite of {\tt Mathematica} packages -comprising a computer algebra toolbox for numerical relativists. {\tt Kranc} -can be used as a ``rapid prototyping'' system for physicists or -mathematicians handling complex systems of partial differential -equations, and through integration into the {\tt Cactus} framework one can -also produce efficient production codes. -We use fourth-order accurate finite differencing for the spacetime -variables and add a fifth-order Kreiss-Oliger dissipation term to -remove high frequency noise. We use a fourth-order Runge-Kutta time -integrator for all evolved variables. \subsubsection{Initial Conditions} Directory: /local_bibtex/ ========================= File [modified]: ott_references.bib Delta lines: +25 -0 =================================================================== --- local_bibtex/ott_references.bib 2011-04-27 17:58:18 UTC (rev 66) +++ local_bibtex/ott_references.bib 2011-04-27 18:24:03 UTC (rev 67) @@ -115,3 +115,28 @@ adsurl = {http://adsabs.harvard.edu/abs/2010RvMP...82.3069C}, } +@Article{Muller:2009jx, + author = {M{\"u}ller, Doreen and Br{\"u}gmann, Bernd}, + title = {Toward a dynamical shift condition for unequal mass + black hole binary simulations}, + journal = {Class. Quantum Grav.}, + volume = 27, + year = 2010, + pages = 114008, + doi = {10.1088/0264-9381/27/11/114008}, + SLACcitation = "%%CITATION = 0912.3125;%%" +} + +@Article{ES-Schnetter2010a, + status = {refereed}, + author = {Erik Schnetter}, + title = {Time Step Size Limitation Introduced by the {BSSN} + {Gamma} Driver}, + journal = {Class. Quantum Grav.}, + volume = 27, + pages = 167001, + year = 2010, + receiveddate = {2010-03-04}, + doi = {10.1088/0264-9381/27/16/167001}, +} +

1 0

[svn:einsteintoolkit] Paper_EinsteinToolkit_2010/ (Rev. 66)
by cott＠tapir.caltech.edu 27 Apr '11

27 Apr '11

User: cott Date: 2011/04/27 12:58 PM Modified: / ET.tex Log: * more intro text to curvature evolution File Changes: Directory: / ============ File [modified]: ET.tex Delta lines: +12 -6 =================================================================== --- ET.tex 2011-04-27 17:47:39 UTC (rev 65) +++ ET.tex 2011-04-27 17:58:18 UTC (rev 66) @@ -935,10 +935,12 @@ and Centrella et al.~\cite{centrella:10}. GR hydrodynamics has been reviewed by Font~\cite{font:08}. In the following, we assume the reader to be familiar with general relativity. +We assume $G = c = M_\odot = 1$ throughout. The Einstein Toolkit provides code to evolve the Einstein equations \begin{equation} G^{\mu\nu} = 8 \pi T^{\mu\nu}\,, +\label{eq:einstein} \end{equation} in the $3+1$ split, foliating 4D spacetime into sequences of spacelike 3-hypersurfaces (slices) connected by timelike normal vectors. In the @@ -949,12 +951,16 @@ observer\footnote{A normal observer follows a wordline tanget to the unit normal on the 3-hypersurface.} and the shift vector $\beta^i$ that describes how spatial coordinates change from one slice to the -next. +next. $T^{\mu\nu}$ in equation~\ref{eq:einstein} is the stress-energy +tensor, which we choose to be that of an ideal relativistic fluid, +\begin{equation} +T^{\mu\nu} = \rho h u^\mu u^\nu - g^{\mu\nu} P\,\,, +\end{equation} +where $\rho$ is the rest-mass density, $u^\mu$ is the 4-velocity, +$g^{\mu\nu}$ is the 4-metric, and $h = 1 + \epsilon + P/\rho$ is the +relativistic specific enthalpy with $\epsilon$ and $P$ being the +specific internal energy and the pressure, respectively. - -\todo{state Einstein equations, mention conventions, give stress -energy tensor} - \subsubsection{Spacetime Curvature Evolution} The Einstein Toolkit curvature evolution code \codename{McLachlan}~\cite{brown:09} is auto-generated from tensor equations via \codename{Kranc} @@ -967,7 +973,7 @@ The evolved variables are the conformal factor $\Phi$, the conformal -3-metric $\tilde{gamma}_{ij}$, the trace $K$ of the extrinsic curvature, +3-metric $\tilde{\gamma}_{ij}$, the trace $K$ of the extrinsic curvature, the trace free extrinsic curvature $A_{ij}$ and the conformal connection functions $\tilde{\Gamma}^i$. These are defined in terms of the standard ADM 4-metric $g_{ij}$, 3-metric $\gamma{ij}$, and extrinsic

1 0

← Newer
1
2
3
4
5
6
7
...
21
Older →

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits April 2011