aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_ShockTubeM.F90
diff options
context:
space:
mode:
Diffstat (limited to 'src/GRHydro_ShockTubeM.F90')
-rw-r--r--src/GRHydro_ShockTubeM.F90573
1 files changed, 565 insertions, 8 deletions
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90
index 186cf69..72bac3e 100644
--- a/src/GRHydro_ShockTubeM.F90
+++ b/src/GRHydro_ShockTubeM.F90
@@ -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,6 +58,8 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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
@@ -69,18 +76,27 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-03 12:40:25 +0000