diff options
| -rw-r--r-- | par/balsara1d.par | 8 | ||||
| -rw-r--r-- | schedule.ccl | 42 | ||||
| -rw-r--r-- | src/GRHydro_ShockTubeM.F90 | 2433 |
3 files changed, 2370 insertions, 113 deletions
diff --git a/par/balsara1d.par b/par/balsara1d.par index 3d1980f..f688b16 100644 --- a/par/balsara1d.par +++ b/par/balsara1d.par @@ -81,7 +81,7 @@ Cactus::terminate="iteration" #Cactus::cctk_final_time = 0.02 #Cactus::cctk_final_time = 0.4 #cactus::cctk_itlast = 40 -cactus::cctk_itlast = 400 +cactus::cctk_itlast = 800 #cactus::cctk_itlast = 2 IO::out_dir = $parfile @@ -97,12 +97,12 @@ CarpetIOBasic::outInfo_every=1 ##CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::psidc grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec" ##CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel grhydro::dens grhydro::tau grhydro::scon " -#CarpetIOHDF5::out_every = 1 -#CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel HydroBase::w_lorentz grhydro::dens grhydro::divB grhydro::tau grhydro::scon HydroBase::Bvec" +CarpetIOHDF5::out_every = 10 +CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel HydroBase::w_lorentz grhydro::dens grhydro::divB grhydro::tau grhydro::scon HydroBase::Bvec" ##CarpetIOHDF5::out_vars = "HydroBase::rho HydroBase::press HydroBase::eps HydroBase::vel HydroBase::w_lorentz grhydro::psidc grhydro::dens grhydro::divB 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 " -CarpetIOHDF5::out2D_every = 1 +CarpetIOHDF5::out2D_every = 10 CarpetIOHDF5::out2D_xy = "yes" CarpetIOHDF5::out2D_xz = "no" CarpetIOHDF5::out2D_yz = "no" diff --git a/schedule.ccl b/schedule.ccl index 2e32506..2b72ab4 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -43,20 +43,46 @@ if (CCTK_Equals(initial_hydro,"shocktube")) { 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(clean_divergence){ + schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons, GRHydro::psidc + } "Diagonal shock boundary conditions" + } else { + schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons + } "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" + schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection + { + LANG: Fortran + } "2-D Diagonal shock boundary conditions" } +# if(CCTK_Equals(shocktube_type,"diagshock2d")) +# { +# if(clean_divergence){ +# schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection +# { +# LANG: Fortran +# SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons, GRHydro::psidc +# } "2-D Diagonal shock boundary conditions" +# } else { +# schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection +# { +# LANG: Fortran +# SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons +# } "2-D Diagonal shock boundary conditions" +# } +# } + } else { schedule GRHydro_shocktube in HydroBase_Initial { diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90 index 6418baf..fc6a45b 100644 --- a/src/GRHydro_ShockTubeM.F90 +++ b/src/GRHydro_ShockTubeM.F90 @@ -57,7 +57,7 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, nx, ny, nz - CCTK_REAL :: direction, det + CCTK_REAL :: direction, det, DIRECTION_TINY CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, & velzl, velzr, epsl, epsr CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr @@ -531,8 +531,10 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS) 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 + DIRECTION_TINY = 0.0001*CCTK_DELTA_SPACE(1) + + if ( ((change_shock_direction==0).and.(direction .lt. DIRECTION_TINY)).or.& + ((change_shock_direction==1).and.(direction .gt. DIRECTION_TINY)) ) then !!$ Left state @@ -651,10 +653,18 @@ subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew + CCTK_INT :: i, j, k, nx, ny, nz, sten, stenp1, minsum, maxsum + CCTK_INT :: inew, jnew, knew CCTK_INT :: xoff,yoff,zoff,indsum CCTK_REAL :: det + CCTK_INT :: status + character(len=200) :: error_message + integer :: gindex_tau, gindex_dens, gindex_scon, gindex_Bcons + CCTK_INT, parameter :: num_groups = 4 + CCTK_INT :: groups(num_groups) + + sten=GRHydro_stencil stenp1=GRHydro_stencil + 1 nx = cctk_lsh(1) @@ -668,42 +678,879 @@ subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS) 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 +!!$ Loop over the cubical domain 6 faces: + + !!$ 1) xmin: + do k = stenp1, nz-sten + do j = stenp1, ny-sten + do i = 1, sten + + indsum = i+j+k + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + dens(i,j,k) = dens(inew,jnew,knew) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + 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) + end do + end do + end do + + !!$ 2) xmax: + do k = stenp1, nz-sten + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + indsum = i+j+k + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners else - yoff=0 + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + dens(i,j,k) = dens(inew,jnew,knew) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) 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 + end do + end do + end do + + !!$ 3) ymin: + do k = stenp1, nz-sten + do j = 1, sten + do i = stenp1, nx-sten + + indsum = i+j+k + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + dens(i,j,k) = dens(inew,jnew,knew) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif - indsum = i+j+k + end do + end do + end do + + !!$ 4) ymax: + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + indsum = i+j+k + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + dens(i,j,k) = dens(inew,jnew,knew) + end if - 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. + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + +!!$ the following two cases are different for a 3d diagonal than a 2d diagonal + + !!$ 5) zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = stenp1, nx-sten + + indsum = i+j+k + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + dens(i,j,k) = dens(inew,jnew,knew) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 6) zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = stenp1, nx-sten + + indsum = i+j+k + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + dens(i,j,k) = dens(inew,jnew,knew) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + +!!$ Loop over the cubical domain 12 edges: + + !!$ 1) xmin,ymin: + do k = stenp1, nz-sten + do j = 1, sten + do i = 1, sten + + indsum = i+j+k +!!$ indsum > maxsum impossible! + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + end if + +!!$ Don't poison edges! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 2) xmin,ymax: + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = 1, sten + + indsum = i+j+k +!!$ indsum > maxsum and indsum < minsum both impossible! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 3) xmax,ymin: + do k = stenp1, nz-sten + do j = 1, sten + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum > maxsum and indsum < minsum both impossible! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 4) xmax,ymax: + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum < minsum impossible! + if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + end if + +!!$ Don't poison edges! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 5) ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = stenp1, nx-sten + + indsum = i+j+k +!!$ indsum > maxsum impossible! + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + end if + +!!$ Don't poison edges! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 6) ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = stenp1, nx-sten + + indsum = i+j+k +!!$ indsum > maxsum and indsum < minsum both impossible! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 7) ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + indsum = i+j+k +!!$ indsum > maxsum and indsum < minsum both impossible! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 8) ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + indsum = i+j+k +!!$ indsum < minsum impossible! + if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + end if + +!!$ Don't poison edges! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 9) xmin,zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = 1, sten + + indsum = i+j+k +!!$ indsum > maxsum impossible! + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + knew = stenp1 + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + end if + +!!$ Don't poison edges! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 10) xmin,zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = 1, sten + + indsum = i+j+k +!!$ indsum > maxsum and indsum < minsum both impossible! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 11) xmax,zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum > maxsum and indsum < minsum both impossible! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 12) xmax,zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum < minsum impossible! + if(indsum.gt.maxsum) then + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + else + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + end if + +!!$ Don't poison edges! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + +!!$ Loop over the cubical domain 8 corners: + + !!$ 1) xmin,ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = 1, sten + +!!$ indsum < minsum guaranteed! + + inew = stenp1 + jnew = stenp1 + knew = stenp1 + +!!$ Don't poison corners! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 2) xmin,ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = 1, sten + + indsum = i+j+k +!!$ indsum < maxsum and indsum > minsum guaranteed! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 3) xmin,ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = 1, sten + + indsum = i+j+k +!!$ indsum < maxsum and indsum > minsum guaranteed! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 4) xmin,ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = 1, sten + + indsum = i+j+k +!!$ indsum < maxsum and indsum > minsum guaranteed! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 5) xmax,ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum < maxsum and indsum > minsum guaranteed! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 6) xmax,ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum < maxsum and indsum > minsum guaranteed! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 7) xmax,ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = nx-sten+1, nx + + indsum = i+j+k +!!$ indsum < maxsum and indsum > minsum guaranteed! + call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + + !!$ 8) xmax,ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = nx-sten+1, nx + +!!$ indsum > maxsum guaranteed! + + inew = nx - stenp1 + 1 + jnew = ny - stenp1 + 1 + knew = nz - stenp1 + 1 + +!!$ Don't poison corners! + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + + end do + end do + end do + +!!$ Synchronize ghost zones (first time): + + call CCTK_GroupIndex(gindex_tau,"GRHydro::tau") + call CCTK_GroupIndex(gindex_dens,"GRHydro::dens") + call CCTK_GroupIndex(gindex_scon,"GRHydro::scon") + call CCTK_GroupIndex(gindex_Bcons,"GRHydro::Bcons") + + groups(1) = gindex_tau + groups(2) = gindex_dens + groups(3) = gindex_scon + groups(4) = gindex_Bcons + + call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups) + if(status.lt.0) then + write(error_message,'(a33,i3,a5,i3,a7)' ) & + 'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + + if(clean_divergence.ne.0) then + call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc") + if(status.lt.0) then + write(error_message,'(a31,i3,a19)' ) & + 'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + endif + + +!!$ Fix the indsum < minsum area in one go? + do k=1,minsum-3 + do j=1,minsum-k-2 + do i=1,minsum-k-j-1 + + indsum=i+j+k + if(dens(i,j,k).lt.0.0d0) then + if(dens(j,i,k).gt.0.0d0) then + inew=j + jnew=i + knew=k + else if(dens(k,j,i).gt.0.0d0) then + inew=k + jnew=j + knew=i + else if(dens(i,k,j).gt.0.0d0) then + inew=i + jnew=k + knew=j + else + inew=stenp1 + jnew=stenp1 + knew=stenp1 + endif - 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,knew) + Bconsy(i,j,k)=Bconsy(inew,jnew,knew) + Bconsz(i,j,k)=Bconsz(inew,jnew,knew) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + end if + + end do + end do + end do + +!!$ Fix the indsum > maxsum area in one go? + do k=nz-3*stenp1+4,nz + do j=ny-3*stenp1+4+(nz-k),ny + do i=nx-3*stenp1+4+(nz-k)+(ny-j),nx + + indsum=i+j+k + xoff=nx-i + yoff=ny-j + zoff=nz-k + if(dens(i,j,k).lt.0.0d0) then + + if(dens(nx-yoff,ny-xoff,k).gt.0.0d0) then + inew=nx-yoff + jnew=ny-xoff + knew=k + else if(dens(nx-zoff,j,nz-xoff).gt.0.0d0) then + inew=nx-zoff + jnew=j + knew=nz-xoff + else if(dens(i,ny-zoff,nz-yoff).gt.0.0d0) then + inew=i + jnew=ny-zoff + knew=nz-xoff + else + inew=nx - stenp1 + 1 + jnew=ny - stenp1 + 1 + knew=nz - stenp1 + 1 + endif dens(i,j,k) = dens(inew,jnew,knew) sx(i,j,k) = sx(inew,jnew,knew) @@ -713,16 +1560,214 @@ subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS) Bconsx(i,j,k)=Bconsx(inew,jnew,knew) Bconsy(i,j,k)=Bconsy(inew,jnew,knew) Bconsz(i,j,k)=Bconsz(inew,jnew,knew) - if(clean_divergence.ne.0) then - psidc(i,j,k)=psidc(inew,jnew,knew) - endif - - endif + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,knew) + endif + end if - enddo - enddo - enddo + end do + end do + end do + +!!$ Fix the ymin face in the region indsum < minsum +!!$ do k = stenp1, nz-sten +!!$ do j = 1, sten +!!$ do i = stenp1, nx-sten +!!$ +!!$ indsum = i+j +!!$ if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then +!!$ inew = j +!!$ jnew = i +!!$ if(dens(inew,jnew,k).lt.0.0d0) then +!!$ inew = stenp1 +!!$ jnew = stenp1 +!!$ endif +!!$ +!!$ 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) +!!$ Bconsx(i,j,k)=Bconsx(inew,jnew,k) +!!$ Bconsy(i,j,k)=Bconsy(inew,jnew,k) +!!$ Bconsz(i,j,k)=Bconsz(inew,jnew,k) +!!$ if(clean_divergence.ne.0) then +!!$ psidc(i,j,k)=psidc(inew,jnew,k) +!!$ endif +!!$ end if +!!$ +!!$ end do +!!$ end do +!!$ end do + +!!$ Fix the xmin face in the region indsum < minsum +!!$ do k = stenp1, nz-sten +!!$ do j = stenp1, ny-sten +!!$ do i = 1, sten +!!$ +!!$ indsum = i+j +!!$ if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then +!!$ inew = j +!!$ jnew = i +!!$ if(dens(inew,jnew,k).lt.0.0d0) then +!!$ inew = stenp1 +!!$ jnew = stenp1 +!!$ endif +!!$ +!!$ 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) +!!$ Bconsx(i,j,k)=Bconsx(inew,jnew,k) +!!$ Bconsy(i,j,k)=Bconsy(inew,jnew,k) +!!$ Bconsz(i,j,k)=Bconsz(inew,jnew,k) +!!$ if(clean_divergence.ne.0) then +!!$ psidc(i,j,k)=psidc(inew,jnew,k) +!!$ endif +!!$ end if +!!$ +!!$ +!!$ end do +!!$ end do +!!$ end do + +!!$ Fix the ymax face in the region indsum > maxsum +!!$ do k = stenp1, nz-sten +!!$ do j = ny-sten+1, ny +!!$ do i = stenp1, nx-sten +!!$ +!!$ indsum = i+j +!!$ if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then +!!$ inew = j-ny+nx +!!$ jnew = i-nx+ny +!!$ if(dens(inew,jnew,k).lt.0.0d0) then +!!$ inew = nx - stenp1 +!!$ jnew = ny - stenp1 +!!$ endif +!!$ +!!$ 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) +!!$ Bconsx(i,j,k)=Bconsx(inew,jnew,k) +!!$ Bconsy(i,j,k)=Bconsy(inew,jnew,k) +!!$ Bconsz(i,j,k)=Bconsz(inew,jnew,k) +!!$ if(clean_divergence.ne.0) then +!!$ psidc(i,j,k)=psidc(inew,jnew,k) +!!$ endif +!!$ endif +!!$ +!!$ end do +!!$ end do +!!$ end do + +!!$ Fix the xmax face in the region indsum > maxsum +!!$ do k = stenp1, nz-sten +!!$ do j = stenp1, ny-sten +!!$ do i = nx-sten+1, nx +!!$ +!!$ indsum = i+j +!!$ if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then +!!$ inew = j-ny+nx +!!$ jnew = i-nx+ny +!!$ if(dens(inew,jnew,k).lt.0.0d0) then +!!$ inew = nx - stenp1 +!!$ jnew = ny - stenp1 +!!$ endif +!!$ +!!$ 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) +!!$ Bconsx(i,j,k)=Bconsx(inew,jnew,k) +!!$ Bconsy(i,j,k)=Bconsy(inew,jnew,k) +!!$ Bconsz(i,j,k)=Bconsz(inew,jnew,k) +!!$ if(clean_divergence.ne.0) then +!!$ psidc(i,j,k)=psidc(inew,jnew,k) +!!$ endif +!!$ endif +!!$ +!!$ end do +!!$ end do +!!$ end do + +!!$ Synchronize ghost zones (second time): + + call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups) + if(status.lt.0) then + write(error_message,'(a33,i3,a5,i3,a7)' ) & + 'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + + if(clean_divergence.ne.0) then + call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc") + if(status.lt.0) then + write(error_message,'(a31,i3,a19)' ) & + 'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + endif + +!!$ 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) +!!$ Bconsx(i,j,k)=Bconsx(inew,jnew,knew) +!!$ Bconsy(i,j,k)=Bconsy(inew,jnew,knew) +!!$ Bconsz(i,j,k)=Bconsz(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) @@ -733,10 +1778,16 @@ subroutine GRHydro_Diagshock2D_BoundaryM(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_INT :: i, j, k, kc, nx, ny, nz, sten, stenp1, minsum, maxsum, inew, jnew, knew + CCTK_INT :: xoff,yoff,zoff,indsum,numoff CCTK_REAL :: det + CCTK_INT :: status + character(len=200) :: error_message + integer :: gindex_tau, gindex_dens, gindex_scon, gindex_Bcons + CCTK_INT, parameter :: num_groups = 4 + CCTK_INT :: groups(num_groups) + sten=GRHydro_stencil stenp1=GRHydro_stencil + 1 nx = cctk_lsh(1) @@ -744,80 +1795,1260 @@ subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS) nz = cctk_lsh(3) minsum = 2*stenp1 - maxsum = nx+ny-2*stenp1 + 2 + maxsum = nx+ny-2*sten 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) +!!$ Loop over the cubical domain 6 faces: + + !!$ 1) xmin: + do k = stenp1, nz-sten + do j = stenp1, ny-sten + do i = 1, sten + + xoff=i-stenp1 + yoff=0 + zoff=0 + indsum = i+j + inew = i - xoff + jnew = j + xoff + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - sten + jnew = ny - sten + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners else - yoff=0 + dens(i,j,k) = dens(inew,jnew,k) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) 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 - + end do + end do + end do + + !!$ 2) xmax: + do k = stenp1, nz-sten + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + xoff=i-(nx-stenp1+1) + yoff=0 + zoff=0 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 = i - xoff + jnew = j + xoff + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - sten + jnew = ny - sten + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + dens(i,j,k) = dens(inew,jnew,k) + end if - inew=indsum/2 - jnew=indsum-inew + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do + + !!$ 3) ymin: + do k = stenp1, nz-sten + do j = 1, sten + do i = stenp1, nx-sten + + xoff=0 + yoff=j-stenp1 + zoff=0 + indsum = i+j + inew = i + yoff + jnew = j - yoff + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - sten + jnew = ny - sten + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + dens(i,j,k) = dens(inew,jnew,k) + end if - 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) - Bconsx(i,j,k)=Bconsx(inew,jnew,k) - Bconsy(i,j,k)=Bconsy(inew,jnew,k) - Bconsz(i,j,k)=Bconsz(inew,jnew,k) - if(clean_divergence.ne.0) then - psidc(i,j,k)=psidc(inew,jnew,k) - endif + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do - else if( zoff.ne.0) then + !!$ 4) ymax: + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + xoff=0 + yoff=j-(ny-stenp1+1) + zoff=0 + indsum = i+j + inew = i + yoff + jnew = j - yoff + if(indsum.lt.minsum) then + inew = stenp1 + jnew = stenp1 + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else if(indsum.gt.maxsum) then + inew = nx - sten + jnew = ny - sten + dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners + else + dens(i,j,k) = dens(inew,jnew,k) + end if + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do + + !!$ 5) zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = stenp1, nx-sten - kc = nz/2 + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif - 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) - Bconsx(i,j,k)=Bconsx(i,j,kc) - Bconsy(i,j,k)=Bconsy(i,j,kc) - Bconsz(i,j,k)=Bconsz(i,j,kc) - if(clean_divergence.ne.0) then - psidc(i,j,k)=psidc(i,j,kc) - endif + end do + end do + end do + + !!$ 6) zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + +!!$ Loop over the cubical domain 12 edges: + + !!$ 1) xmin,ymin: + do k = stenp1, nz-sten + do j = 1, sten + do i = 1, sten + + !!$ We already know that for this edge indsum.lt.minsum. + inew = stenp1 + jnew = stenp1 + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do + + !!$ 2) xmin,ymax: + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = 1, sten + + xoff=i-stenp1 + yoff=j-(ny-stenp1+1) + zoff=0 + + numoff = max(abs(xoff),abs(yoff)) + if( abs(xoff).eq.numoff) then + inew = i + numoff + jnew = j - numoff + else + jnew = j - numoff + inew = i + numoff + endif + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do + + !!$ 3) xmax,ymin: + do k = stenp1, nz-sten + do j = 1, sten + do i = nx-sten+1, nx + + xoff=i-(nx-stenp1+1) + yoff=j-stenp1 + zoff=0 + numoff = max(abs(xoff),abs(yoff)) + if( abs(xoff).eq.numoff) then + inew = i - numoff + jnew = j + numoff + else + jnew = j + numoff + inew = i - numoff + endif + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do + + !!$ 4) xmax,ymax: + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = nx-sten+1, nx + + !!$ We already know that for this edge indsum.gt.maxsum + inew = nx - sten + jnew = ny - sten + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + + end do + end do + end do + + !!$ 5) ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 6) ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 7) ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 8) ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 9) xmin,zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 10) xmin,zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 11) xmax,zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 12) xmax,zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + +!!$ Loop over the cubical domain 8 corners: + + !!$ 1) xmin,ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 2) xmin,ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 3) xmin,ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 4) xmin,ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 5) xmax,ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 6) xmax,ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 7) xmax,ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 8) xmax,ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + +!!$ Synchronize ghost zones (first time): + + call CCTK_GroupIndex(gindex_tau,"GRHydro::tau") + call CCTK_GroupIndex(gindex_dens,"GRHydro::dens") + call CCTK_GroupIndex(gindex_scon,"GRHydro::scon") + call CCTK_GroupIndex(gindex_Bcons,"GRHydro::Bcons") + + groups(1) = gindex_tau + groups(2) = gindex_dens + groups(3) = gindex_scon + groups(4) = gindex_Bcons + + call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups) + if(status.lt.0) then + write(error_message,'(a33,i3,a5,i3,a7)' ) & + 'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + + if(clean_divergence.ne.0) then + call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc") + if(status.lt.0) then + write(error_message,'(a31,i3,a19)' ) & + 'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + endif + +!!$ Fix the ymin face in the region indsum < minsum + do k = stenp1, nz-sten + do j = 1, sten + do i = stenp1, nx-sten + + indsum = i+j + if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then + inew = j + jnew = i + if(dens(inew,jnew,k).lt.0.0d0) then + inew = stenp1 + jnew = stenp1 + endif + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + end if + + end do + end do + end do + +!!$ Fix the xmin face in the region indsum < minsum + do k = stenp1, nz-sten + do j = stenp1, ny-sten + do i = 1, sten + + indsum = i+j + if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then + inew = j + jnew = i + if(dens(inew,jnew,k).lt.0.0d0) then + inew = stenp1 + jnew = stenp1 + endif + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + end if + + + end do + end do + end do + +!!$ Fix the ymax face in the region indsum > maxsum + do k = stenp1, nz-sten + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + indsum = i+j + if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then + inew = j-ny+nx + jnew = i-nx+ny + if(dens(inew,jnew,k).lt.0.0d0) then + inew = nx - sten + jnew = ny - sten + endif + + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + endif + + end do + end do + end do + +!!$ Fix the xmax face in the region indsum > maxsum + do k = stenp1, nz-sten + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + indsum = i+j + if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then + inew = j-ny+nx + jnew = i-nx+ny + if(dens(inew,jnew,k).lt.0.0d0) then + inew = nx - sten + jnew = ny - sten + endif + 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) + Bconsx(i,j,k)=Bconsx(inew,jnew,k) + Bconsy(i,j,k)=Bconsy(inew,jnew,k) + Bconsz(i,j,k)=Bconsz(inew,jnew,k) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(inew,jnew,k) + endif + endif + + end do + end do + end do + +!!$ Fix the upper and lower (5-12 below) domain edges since +!!$ they still contain poison got from the faces in the initial +!!$ step: + + !!$ 5) ymin,zmin: + do k = 1, sten + do j = 1, sten + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 6) ymin,zmax: + do k = nz-sten+1, nz + do j = 1, sten + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 7) ymax,zmin: + do k = 1, sten + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 8) ymax,zmax: + do k = nz-sten+1, nz + do j = ny-sten+1, ny + do i = stenp1, nx-sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 9) xmin,zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 10) xmin,zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = 1, sten + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 11) xmax,zmin: + do k = 1, sten + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) + endif + + end do + end do + end do + + !!$ 12) xmax,zmax: + do k = nz-sten+1, nz + do j = stenp1, ny-sten + do i = nx-sten+1, nx + + 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) + Bconsx(i,j,k)=Bconsx(i,j,kc) + Bconsy(i,j,k)=Bconsy(i,j,kc) + Bconsz(i,j,k)=Bconsz(i,j,kc) + if(clean_divergence.ne.0) then + psidc(i,j,k)=psidc(i,j,kc) endif + + end do + end do + end do + +!!$ Synchronize ghost zones (second time): + + call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups) + if(status.lt.0) then + write(error_message,'(a33,i3,a5,i3,a7)' ) & + 'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + + if(clean_divergence.ne.0) then + call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc") + if(status.lt.0) then + write(error_message,'(a31,i3,a19)' ) & + 'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc' + call CCTK_WARN(CCTK_WARN_ABORT,error_message) + endif + endif - enddo - enddo - enddo +!!$ 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) then +!!$ +!!$ if(indsum.ge.minsum.and.indsum.le.maxsum) then +!!$ +!!$ numoff = max(abs(xoff),abs(yoff)) +!!$ if( abs(xoff).eq.numoff) then +!!$ if(xoff.gt.0) then +!!$ inew = i - numoff +!!$ jnew = j + numoff +!!$ else +!!$ inew = i + numoff +!!$ jnew = j - numoff +!!$ endif +!!$ else +!!$ if(yoff.gt.0) then +!!$ jnew = j - numoff +!!$ inew = i + numoff +!!$ else +!!$ jnew = j + numoff +!!$ inew = i - numoff +!!$ endif +!!$ endif +!!$!!$ Note that the following two conditions helps but +!!$!!$ don't solve the problem when running in several +!!$!!$ processors +!!$ else if(indsum.lt.minsum) then +!!$ inew = stenp1 +!!$ jnew = stenp1 +!!$ +!!$ else if(indsum.gt.maxsum) then +!!$ inew = nx - stenp1 +!!$ jnew = ny - stenp1 +!!$ end if +!!$ +!!$ 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) +!!$ Bconsx(i,j,k)=Bconsx(inew,jnew,k) +!!$ Bconsy(i,j,k)=Bconsy(inew,jnew,k) +!!$ Bconsz(i,j,k)=Bconsz(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) +!!$ Bconsx(i,j,k)=Bconsx(i,j,kc) +!!$ Bconsy(i,j,k)=Bconsy(i,j,kc) +!!$ Bconsz(i,j,k)=Bconsz(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 + +subroutine find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew) + +!!$ Find the nearest mapped point in the interior, assuming it exists + + implicit none + CCTK_INT :: i,j,k,n,nx,ny,nz,stenp1 + CCTK_INT :: xoff,yoff,zoff,indsum,minsum,maxsum + CCTK_INT :: inew,jnew,knew,dspace,newsum + + minsum = 3*stenp1 + maxsum = nx+ny+nz - 3*(stenp1-1) + + indsum = i+j+k + + if(indsum.lt.minsum.or.indsum.gt.maxsum)write(6,*)'HUGE ERROR - ILLEGAL INDEX SUM!!!!!!' + + 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 + + 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 + + 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 + + if(xoff.eq.0.and.yoff.eq.0.and.zoff.eq.0) then + write(6,*)'Why did you try to map an interior point?' + inew=i + jnew=j + knew=k + return + else + +!!$ We have something to actually do! + + inew = i-xoff + jnew = j-yoff + knew = k-zoff + dspace = xoff + yoff + zoff + + do n=1,abs(dspace) + + if(dspace.lt.0) then + +!!$ dspace negative; point sum corrected upward +!!$ need to move other coords down + + if(inew.gt.stenp1) then + inew = inew - 1 + else if(jnew.gt.stenp1) then + jnew = jnew - 1 + else if (knew.gt.stenp1) then + knew = knew - 1 + else + write(6,*)'Big error!!!!' + endif + + else if (dspace.gt.0) then + +!!$ dspace positive; point sum corrected downward +!!$ need to move other coords up + + if(inew.lt.nx-stenp1+1) then + inew=inew+1 + else if(jnew.lt.ny-stenp1+1) then + jnew=jnew+1 + else if(knew.lt.nz-stenp1+1) then + knew=knew+1 + else + write(6,*)'Big error2!!!!' + endif + + end if + + end do + + newsum = inew+jnew+knew + + if(inew.lt.stenp1.or.inew.gt.(nx-(stenp1-1)).or. & + jnew.lt.stenp1.or.jnew.gt.(ny-(stenp1-1)).or. & + knew.lt.stenp1.or.knew.gt.(nz-(stenp1-1)).or. & + newsum.ne.indsum)write(6,*)'Big error 3!!!!' + + end if + + return + +end subroutine find_mapped_point + |