diff options
Diffstat (limited to 'src/GRHydro_ShockTubeM.F90')
| -rw-r--r-- | src/GRHydro_ShockTubeM.F90 | 101 |
1 files changed, 93 insertions, 8 deletions
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90 index 1169fef..96b0e7c 100644 --- a/src/GRHydro_ShockTubeM.F90 +++ b/src/GRHydro_ShockTubeM.F90 @@ -54,6 +54,13 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) velzl, velzr, epsl, epsr CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr + bvcxl = Bx_init + bvcyl = By_init + bvczl = Bz_init + bvcxr = Bx_init + bvcyr = By_init + bvczr = Bz_init + nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) @@ -110,17 +117,95 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) velzr = 0.d0 epsl = 1500.d0 epsr = 1.5d-2 - else + else if (CCTK_EQUALS(shock_case,"Balsara1")) then + rhol = 1.d0 + rhor = 0.125d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + bvcxl=0.5d0 + bvcxr=0.5d0 + bvcyl=1.d0 + bvcyr=-1.d0 + bvczl=0.d0 + bvczr=0.d0 + epsl = 1.0/rhol + epsr = 0.1/rhor + else if (CCTK_EQUALS(shock_case,"Balsara2")) then + rhol = 1.d0 + rhor = 1.d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + bvcxl=5.0d0 + bvcxr=5.0d0 + bvcyl=6.d0 + bvcyr=0.7d0 + bvczl=6.d0 + bvczr=0.7d0 + epsl = 1.5d0*30.0d0/rhol + epsr = 1.5d0*1.0d0/rhor + else if (CCTK_EQUALS(shock_case,"Balsara3")) then + rhol = 1.d0 + rhor = 1.d0 + velxl = 0.d0 + velxr = 0.d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + bvcxl=10.0d0 + bvcxr=10.0d0 + bvcyl=7.d0 + bvcyr=0.7d0 + bvczl=7.d0 + bvczr=0.7d0 + epsl = 1.5d0*1000.0d0/rhol + epsr = 1.5d0*0.1d0/rhor + else if (CCTK_EQUALS(shock_case,"Balsara4")) then + rhol = 1.d0 + rhor = 1.d0 + velxl = 0.999d0 + velxr = -0.999d0 + velyl = 0.d0 + velyr = 0.d0 + velzl = 0.d0 + velzr = 0.d0 + bvcxl=10.0d0 + bvcxr=10.0d0 + bvcyl=7.d0 + bvcyr=-7.d0 + bvczl=7.d0 + bvczr=-7.d0 + epsl = 1.5d0*0.1d0/rhol + epsr = 1.5d0*0.1d0/rhor + else if (CCTK_EQUALS(shock_case,"Balsara5")) then + rhol = 1.08d0 + rhor = 1.d0 + velxl = 0.4d0 + velxr = -0.45d0 + velyl = 0.3d0 + velyr = -0.2d0 + velzl = 0.2d0 + velzr = 0.2d0 + bvcxl=2.0d0 + bvcxr=2.0d0 + bvcyl=0.3d0 + bvcyr=-0.7d0 + bvczl=0.3d0 + bvczr=0.5d0 + epsl = 1.5d0*0.95d0/rhol + epsr = 1.5d0*1.0d0/rhor + else call CCTK_WARN(0,"Shock case not recognized") end if - bvcxl = Bx_init - bvcyl = By_init - bvczl = Bz_init - bvcxr = Bx_init - bvcyr = By_init - bvczr = Bz_init - if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then rho(i,j,k) = rhol |