First try of a constraint transport scheme for unigrid

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@296 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

6 files changed, 166 insertions, 3 deletions

diff --git a/interface.ccl b/interface.ccl
index 2b58ca5..c69e3b2 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -347,6 +347,8 @@ real scon[3] type = GF Timelevels = 3 tags='ProlongationParameter="HydroBase::pr
 
 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"'
diff --git a/param.ccl b/param.ccl
index d750946..a4f037c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -537,6 +537,10 @@ boolean track_divB "Track the magnetic field constraint violations"
 {
 } "no"
 
+boolean transport_constraints "Use constraint transport for magnetic field"
+{
+} "no"
+
 boolean calculate_bcom  "Calculate the comoving contravariant magnetic 4-vector b^a?"
 {
 } "no"
diff --git a/schedule.ccl b/schedule.ccl
index 05d031b..ad870ce 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -107,6 +107,11 @@ if (evolve_tracer)
   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
diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90
index fffb42d..cc7b136 100644
--- a/src/GRHydro_CalcUpdate.F90
+++ b/src/GRHydro_CalcUpdate.F90
@@ -48,9 +48,9 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
     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 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
                        (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 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
   
   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 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
     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
 
diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90
index 6758c90..7b04b8c 100644
--- a/src/GRHydro_Source.F90
+++ b/src/GRHydro_Source.F90
@@ -93,6 +93,11 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
   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 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
 
   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
 
 
diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90
index 73334e3..1456500 100644
--- a/src/GRHydro_SourceM.F90
+++ b/src/GRHydro_SourceM.F90
@@ -127,6 +127,13 @@ subroutine SourceTermsM(CCTK_ARGUMENTS)
   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