second iteration of constraint transport

  * fix some indices
  * move poison loop to proper spot
  * compute conserved divergence in track_divB
  * add Maple worksheet to check constraint transport indices


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

1 files changed, 61 insertions, 47 deletions

diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90
index cc7b136..e53f4b8 100644
--- a/src/GRHydro_CalcUpdate.F90
+++ b/src/GRHydro_CalcUpdate.F90
@@ -39,7 +39,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
 
   CCTK_INT :: i,j,k,itracer
-  CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r
+  CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r, alp_tmp
 
   idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction)
 
@@ -73,39 +73,58 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
                  (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - &
                  alp_r * tauflux(i,j,k)) * idx
             if(evolve_mhd.ne.0) then
-               Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + &
-                    (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - &
-                    alp_r * Bconsxflux(i,j,k)) * idx
-               Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + &
-                    (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - &
-                    alp_r * Bconsyflux(i,j,k)) * idx
-               Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + &
-                    (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - &
-                    alp_r * Bconszflux(i,j,k)) * idx
+               if(transport_constraints.ne.0) then
+                  ! we have to first compute all components of v\crossB = E and
+                  ! combine them in the last substep into Bconshs
+                  ! 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
+                     alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1))
+                     Evec(i,j,k,2) = Evec(i,j,k,2) + 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i  ,j  ,k+1))
+                     alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset))
+                     Evec(i,j,k,3) = Evec(i,j,k,3) - 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i  ,j+1,k  ))
+                  elseif(flux_direction.eq.2) then
+                     alp_tmp = 0.5d0 * (alp(i,j,k+1) + alp(i+xoffset,j+yoffset,k+zoffset+1))
+                     Evec(i,j,k,1) = Evec(i,j,k,1) - 0.25d0 * (alp_r*Bconszflux(i,j,k) + alp_tmp*Bconszflux(i  ,j  ,k+1))
+                     alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset))
+                     Evec(i,j,k,3) = Evec(i,j,k,3) + 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j  ,k  ))
+                  elseif(flux_direction.eq.3) then
+                     alp_tmp = 0.5d0 * (alp(i,j+1,k) + alp(i+xoffset,j+yoffset+1,k+zoffset))
+                     Evec(i,j,k,1) = Evec(i,j,k,1) + 0.25d0 * (alp_r*Bconsyflux(i,j,k) + alp_tmp*Bconsyflux(i  ,j+1,k  ))
+                     alp_tmp = 0.5d0 * (alp(i+1,j,k) + alp(i+xoffset+1,j+yoffset,k+zoffset))
+                     Evec(i,j,k,2) = Evec(i,j,k,2) - 0.25d0 * (alp_r*Bconsxflux(i,j,k) + alp_tmp*Bconsxflux(i+1,j  ,k  ))
+                  end if
+               else
+                  Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + &
+                       (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - &
+                       alp_r * Bconsxflux(i,j,k)) * idx
+                  Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + &
+                       (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - &
+                       alp_r * Bconsyflux(i,j,k)) * idx
+                  Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + &
+                       (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - &
+                       alp_r * Bconszflux(i,j,k)) * idx
+               endif
                if(clean_divergence.ne.0) then
                   psidcrhs(i,j,k) = psidcrhs(i,j,k) + &
                        (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))
-                 Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + &
-                          Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction))
-                 divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx 
+                 if(transport_constraints.ne.0) then
+                   ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27)
+                   ! FIXME: this uses the Bcons before the current MoL substep update, should be in MoL_PseudoEvolution instead
+                   divB(i,j,k) = divB(i,j,k) + &
+                     0.25d0*(Bcons(i+xoffset,j+yoffset  ,k+zoffset,flux_direction)-Bcons(i          ,j          ,k          ,flux_direction)+ &
+                             Bcons(i+1      ,j+1-zoffset,k+zoffset,flux_direction)-Bcons(i+1-xoffset,j+xoffset  ,k          ,flux_direction)+ &
+                             Bcons(i+xoffset,j+1-xoffset,k+1      ,flux_direction)-Bcons(i          ,j+zoffset  ,k+1-zoffset,flux_direction)+ &
+                             Bcons(i+1      ,j+1        ,k+1      ,flux_direction)-Bcons(i+1-xoffset,j+1-yoffset,k+1-zoffset,flux_direction))*idx
+                 else
+                   Bvec_l = 0.5d0 * (Bvec(i,j,k,flux_direction) + &
+                            Bvec(i-xoffset,j-yoffset,k-zoffset,flux_direction))
+                   Bvec_r = 0.5d0 * (Bvec(i,j,k,flux_direction) + &
+                            Bvec(i+xoffset,j+yoffset,k+zoffset,flux_direction))
+                   divB(i,j,k) = divB(i,j,k) + ( alp_l * Bvec_l - alp_r * Bvec_r ) * idx 
+                endif
                endif
 
             endif
@@ -272,32 +291,27 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
   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
+    ! FIXME: 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)
+          Bconsrhs(i,j,k,1) = - 0.5d0 * (Evec(i-1,j  ,k-1,2)-Evec(i-1,j  ,k  ,2)) / CCTK_DELTA_SPACE(3) &
+                              - 0.5d0 * (Evec(i-1,j  ,k  ,3)-Evec(i-1,j-1,k  ,3)) / CCTK_DELTA_SPACE(2) &
+                              - 0.5d0 * (Evec(i  ,j  ,k-1,2)-Evec(i  ,j  ,k  ,2)) / CCTK_DELTA_SPACE(3) &
+                              - 0.5d0 * (Evec(i  ,j  ,k  ,3)-Evec(i  ,j-1,k  ,3)) / CCTK_DELTA_SPACE(2)
+          Bconsrhs(i,j,k,2) = - 0.5d0 * (Evec(i-1,j-1,k  ,3)-Evec(i  ,j-1,k  ,3)) / CCTK_DELTA_SPACE(1) &
+                              - 0.5d0 * (Evec(i  ,j-1,k  ,1)-Evec(i  ,j-1,k-1,1)) / CCTK_DELTA_SPACE(3) &
+                              - 0.5d0 * (Evec(i-1,j  ,k  ,3)-Evec(i  ,j  ,k  ,3)) / CCTK_DELTA_SPACE(1) &
+                              - 0.5d0 * (Evec(i  ,j  ,k  ,1)-Evec(i  ,j  ,k-1,1)) / CCTK_DELTA_SPACE(3)
+          Bconsrhs(i,j,k,3) = - 0.5d0 * (Evec(i  ,j-1,k-1,1)-Evec(i  ,j  ,k-1,1)) / CCTK_DELTA_SPACE(2) &
+                              - 0.5d0 * (Evec(i  ,j  ,k-1,2)-Evec(i-1,j  ,k-1,2)) / CCTK_DELTA_SPACE(1) &
+                              - 0.5d0 * (Evec(i  ,j-1,k  ,1)-Evec(i  ,j  ,k  ,1)) / CCTK_DELTA_SPACE(2) &
+                              - 0.5d0 * (Evec(i  ,j  ,k  ,2)-Evec(i-1,j  ,k  ,2)) / CCTK_DELTA_SPACE(1)
         enddo
       enddo
     enddo
   end if
-
-          
   
   return