[Commits] [svn:einsteintoolkit] GRHydro_InitData/trunk/ (Rev. 123)
bcmsma at astro.rit.edu
bcmsma at astro.rit.edu
Wed Apr 27 17:25:27 CDT 2011
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
More information about the Commits
mailing list