[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 296)
roland.haas at physics.gatech.edu
roland.haas at physics.gatech.edu
Wed Nov 2 10:04:18 CDT 2011
User: rhaas
Date: 2011/11/02 10:04 AM
Modified:
/trunk/
interface.ccl, param.ccl, schedule.ccl
/trunk/src/
GRHydro_CalcUpdate.F90, GRHydro_Source.F90, GRHydro_SourceM.F90
Log:
First try of a constraint transport scheme for unigrid
File Changes:
Directory: /trunk/src/
======================
File [modified]: GRHydro_CalcUpdate.F90
Delta lines: +48 -3
===================================================================
--- trunk/src/GRHydro_CalcUpdate.F90 2011-11-02 13:44:20 UTC (rev 295)
+++ trunk/src/GRHydro_CalcUpdate.F90 2011-11-02 15:04:18 UTC (rev 296)
@@ -48,9 +48,9 @@
if (use_weighted_fluxes == 0) then
!$OMP PARALLEL DO PRIVATE(i,j,k,itracer,alp_l,alp_r,Bvec_l,Bvec_r)
- do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil
- do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil
- do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil
+ do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! mean we need to compute fluxes for one more point the GRHydro_stencil would indicate
+ do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil
+ do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil
alp_l = 0.5d0 * (alp(i,j,k) + &
alp(i-xoffset,j-yoffset,k-zoffset))
@@ -87,6 +87,19 @@
(alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
alp_r * psidcflux(i,j,k)) * idx
endif
+ if(transport_constraints.ne.0) then
+ ! Evec lives on edges of cell: Evec(i,j,k,1) is at edge i,j+1/2,k+1/2 ie. the lower-front edge of cell (i,j,k)
+ if(flux_direction.eq.1) then
+ Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i,j ,k+1))
+ Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsyflux(i,j+1,k ))
+ elseif(flux_direction.eq.2) then
+ Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * alp_r*(Bconszflux(i,j,k) + Bconszflux(i ,j,k+1))
+ Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j,k ))
+ elseif(flux_direction.eq.3) then
+ Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * alp_r*(Bconsyflux(i,j,k) + Bconsyflux(i ,j+1,k))
+ Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * alp_r*(Bconsxflux(i,j,k) + Bconsxflux(i+1,j ,k))
+ end if
+ end if
if(track_divB.ne.0) then
Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + &
Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction))
@@ -205,6 +218,10 @@
else if (CCTK_EQUALS(method_type, "Flux split FD")) then
+ if (transport_constraints .ne. 0) then
+ call CCTK_WARN(0, "Not supported")
+ end if
+
do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil
do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil
do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil
@@ -253,6 +270,34 @@
enddo
end if
+
+ if (transport_constraints.ne.0 .and. flux_direction.eq.1) then ! HACK: x direction is last
+ ! RH: I think one could wrap all of this into a single do loop and remove the
+ ! Evec storage
+ do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil
+ do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil
+ do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil
+ Bconsrhs(i,j,k,1) = & !Bconsrhs(i,j,k,1) &
+ - 0.5d0 * (Evec(i-1,j-1,k-1,2)-Evec(i-1,j-1,k ,2)) / CCTK_DELTA_SPACE(3) &
+ - 0.5d0 * (Evec(i-1,j ,k-1,3)-Evec(i-1,j-1,k-1,3)) / CCTK_DELTA_SPACE(2) &
+ - 0.5d0 * (Evec(i ,j-1,k-1,2)-Evec(i ,j-1,k ,2)) / CCTK_DELTA_SPACE(3) &
+ - 0.5d0 * (Evec(i ,j ,k-1,3)-Evec(i ,j-1,k-1,3)) / CCTK_DELTA_SPACE(2)
+ Bconsrhs(i,j,k,2) = & !Bconsrhs(i,j,k,2) &
+ - 0.5d0 * (Evec(i-1,j-1,k-1,3)-Evec(i-1,j ,k-1,3)) / CCTK_DELTA_SPACE(1) &
+ - 0.5d0 * (Evec(i ,j-1,k-1,1)-Evec(i-1,j-1,k-1,1)) / CCTK_DELTA_SPACE(3) &
+ - 0.5d0 * (Evec(i-1,j-1,k ,3)-Evec(i-1,j ,k ,3)) / CCTK_DELTA_SPACE(1) &
+ - 0.5d0 * (Evec(i ,j-1,k ,1)-Evec(i-1,j-1,k ,1)) / CCTK_DELTA_SPACE(3)
+ Bconsrhs(i,j,k,3) = & !Bconsrhs(i,j,k,3) &
+ - 0.5d0 * (Evec(i-1,j-1,k-1,1)-Evec(i ,j-1,k-1,1)) / CCTK_DELTA_SPACE(1) &
+ - 0.5d0 * (Evec(i-1,j-1,k ,2)-Evec(i-1,j-1,k-1,2)) / CCTK_DELTA_SPACE(2) &
+ - 0.5d0 * (Evec(i-1,j ,k-1,1)-Evec(i ,j ,k-1,1)) / CCTK_DELTA_SPACE(1) &
+ - 0.5d0 * (Evec(i-1,j ,k ,2)-Evec(i-1,j ,k-1,2)) / CCTK_DELTA_SPACE(2)
+ enddo
+ enddo
+ enddo
+ end if
+
+
return
File [modified]: GRHydro_Source.F90
Delta lines: +100 -0
===================================================================
--- trunk/src/GRHydro_Source.F90 2011-11-02 13:44:20 UTC (rev 295)
+++ trunk/src/GRHydro_Source.F90 2011-11-02 15:04:18 UTC (rev 296)
@@ -93,6 +93,11 @@
CCTK_REAL, DIMENSION(:,:,:), POINTER :: beta1, beta2, beta3
CCTK_REAL, DIMENSION(:,:,:,:), POINTER :: vup
+ ! poison hack
+ CCTK_INT, SAVE :: reflevel = -1
+ CCTK_INT, SAVE :: last_iteration_seen = -1
+ CCTK_INT, SAVE :: mol_substep = -1
+
if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
g11 => gaa
g12 => gab
@@ -449,6 +454,101 @@
deallocate(force_spatial_second_order)
+#if(0) /* poison edges of domain */
+ if(last_iteration_seen .ne. cctk_iteration .or. reflevel .ne. grhydro_reflevel) then
+ last_iteration_seen = cctk_iteration
+ reflevel = grhydro_reflevel
+ mol_substep = 0
+ else
+ mol_substep = mol_substep + 1
+ end if
+ do k = 1, GRHydro_stencil*mol_substep
+ do j = 1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Bcons(i,j,k,1) = -1d100
+ Bcons(i,j,k,2) = -1d100
+ Bcons(i,j,k,3) = -1d100
+ end do
+ end do
+ end do
+ do k = cctk_lsh(3)-GRHydro_stencil*mol_substep+1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Bcons(i,j,k,1) = -1d100
+ Bcons(i,j,k,2) = -1d100
+ Bcons(i,j,k,3) = -1d100
+ end do
+ end do
+ end do
+ do i = 1, GRHydro_stencil*mol_substep
+ do k = 1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Bcons(i,j,k,1) = -1d100
+ Bcons(i,j,k,2) = -1d100
+ Bcons(i,j,k,3) = -1d100
+ end do
+ end do
+ end do
+ do i = cctk_lsh(1)-GRHydro_stencil*mol_substep+1, cctk_lsh(1)
+ do k = 1, cctk_lsh(3)
+ do j = 1, cctk_lsh(2)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Bcons(i,j,k,1) = -1d100
+ Bcons(i,j,k,2) = -1d100
+ Bcons(i,j,k,3) = -1d100
+ end do
+ end do
+ end do
+ do j = 1, GRHydro_stencil*mol_substep
+ do i = 1, cctk_lsh(1)
+ do k = 1, cctk_lsh(3)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Bcons(i,j,k,1) = -1d100
+ Bcons(i,j,k,2) = -1d100
+ Bcons(i,j,k,3) = -1d100
+ end do
+ end do
+ end do
+ do j = cctk_lsh(2)-GRHydro_stencil*mol_substep+1, cctk_lsh(2)
+ do i = 1, cctk_lsh(1)
+ do k = 1, cctk_lsh(3)
+ dens(i,j,k) = -1d100
+ Scon(i,j,k,1) = -1d100
+ Scon(i,j,k,2) = -1d100
+ Scon(i,j,k,3) = -1d100
+ tau(i,j,k) = -1d100
+ Bcons(i,j,k,1) = -1d100
+ Bcons(i,j,k,2) = -1d100
+ Bcons(i,j,k,3) = -1d100
+ end do
+ end do
+ end do
+#endif
+
+
end subroutine SourceTerms
File [modified]: GRHydro_SourceM.F90
Delta lines: +7 -0
===================================================================
--- trunk/src/GRHydro_SourceM.F90 2011-11-02 13:44:20 UTC (rev 295)
+++ trunk/src/GRHydro_SourceM.F90 2011-11-02 15:04:18 UTC (rev 296)
@@ -127,6 +127,13 @@
if (track_divB .ne. 0) then
divB=0.d0
endif
+
+ if (transport_constraints .ne. 0) then
+ Evec = 0.d0
+ Bconsrhsx(i,j,k) = 0.d0
+ Bconsrhsy(i,j,k) = 0.d0
+ Bconsrhsz(i,j,k) = 0.d0
+ endif
!!$ Set up the array for checking the order. We switch to second order
Directory: /trunk/
==================
File [modified]: interface.ccl
Delta lines: +2 -0
===================================================================
--- trunk/interface.ccl 2011-11-02 13:44:20 UTC (rev 295)
+++ trunk/interface.ccl 2011-11-02 15:04:18 UTC (rev 296)
@@ -347,6 +347,8 @@
real Bcons[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorparity=-1 tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "B-field conservative variable"
+real Evec[3] type = GF Timelevels = 1 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="D" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Electric field at edges"
+
real Y_e_con type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar" tensorweight=+1.0 jacobian="inverse_jacobian" interpolator="matter"' "Conserved electron fraction"
real GRHydro_tracers[number_of_tracers] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::prolongation_type" tensortypealias="Scalar"'
File [modified]: param.ccl
Delta lines: +4 -0
===================================================================
--- trunk/param.ccl 2011-11-02 13:44:20 UTC (rev 295)
+++ trunk/param.ccl 2011-11-02 15:04:18 UTC (rev 296)
@@ -537,6 +537,10 @@
{
} "no"
+boolean transport_constraints "Use constraint transport for magnetic field"
+{
+} "no"
+
boolean calculate_bcom "Calculate the comoving contravariant magnetic 4-vector b^a?"
{
} "no"
File [modified]: schedule.ccl
Delta lines: +5 -0
===================================================================
--- trunk/schedule.ccl 2011-11-02 13:44:20 UTC (rev 295)
+++ trunk/schedule.ccl 2011-11-02 15:04:18 UTC (rev 296)
@@ -107,6 +107,11 @@
STORAGE:GRHydro_tracer_rhs
}
+if (transport_constraints)
+{
+ STORAGE:Evec
+}
+
# leave this here for compatibility with older code, should not be used anymore
# and go out at some point
schedule group GRHydro_Initial IN HydroBase_Initial BEFORE SetTmunu
More information about the Commits
mailing list