RIT MHD development:

1) Fix a bug in the divergence cleaning implementation:
a psidc term in the induction equation was implemented 
as a source term where it was supposed to be coded as 
part of the flux calculation.

2) Introduce divB as a diagnostic grid function.

3) Remove an old file: GRHydro_CalcUpdateM.F90



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

9 files changed, 133 insertions, 312 deletions

diff --git a/interface.ccl b/interface.ccl
index 2b4439c..de24a28 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -332,6 +332,8 @@ real Bvecrhs[3] type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="n
 
 real psidcrhs type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Update term for psidc"
 
+real divB type = GF Timelevels = 1 tags='Prolongation="None" checkpoint="no"' "Magnetic field constraint"
+
 ##################################################
 ### These variables are only protected so that ###
 ### the tests in init_data work. Should fix.   ###
diff --git a/param.ccl b/param.ccl
index ac01a95..56946aa 100644
--- a/param.ccl
+++ b/param.ccl
@@ -492,7 +492,7 @@ CCTK_REAL GRHydro_lorentz_overshoot_cutoff "Set the Lorentz factor to this value
   0:* :: "Some big value"
 } 1.e100
 
-boolean clean_divergence "Should we advect tracers?"
+boolean clean_divergence "Use hyperbolic divergence cleaning"
 {
 } "no"
 
@@ -506,6 +506,10 @@ CCTK_REAL cp_dc "The c_p parameter for divergence cleaning"
   0:* :: "Any value, but one to 12 is preferred"
 } 1.0
 
+boolean track_divB "Track the magnetic field constraint violations"
+{
+} "no"
+
 ############### temporary parameters to be removed once connected issues are fixed.
 
 boolean con2prim_oct_hack "Disregard c2p failures in oct/rotsym90 boundary cells with xyz < 0" STEERABLE=ALWAYS
diff --git a/schedule.ccl b/schedule.ccl
index 284c773..d9e8552 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -68,10 +68,17 @@ STORAGE:GRHydro_reflevel
 STORAGE:densrhs
 STORAGE:taurhs
 STORAGE:srhs
-STORAGE:Bvecrhs
-if (clean_divergence)
+if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
 {
-  STORAGE:psidcrhs
+  STORAGE:Bvecrhs
+  if (clean_divergence)
+  {
+    STORAGE:psidcrhs
+  }
+  if (track_divB)
+  {
+    STORAGE:divB
+  }
 }
 STORAGE:GRHydro_eos_scalars
 STORAGE:GRHydro_minima
@@ -97,12 +104,22 @@ schedule group GRHydro_Initial IN HydroBase_Initial BEFORE SetTmunu
 {
 } "GRHydro initial data group"
 
-if(CCTK_Equals(Bvec_evolution_method,"GRHydro") && clean_divergence)
+if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
 {
-  schedule GRHydro_InitDivergenceClean IN HydroBase_Initial
+  if (clean_divergence)
   {
-    LANG: Fortran
-  } "Set psi for divergence cleaning initially to zero"
+    schedule GRHydro_InitDivergenceClean IN HydroBase_Initial
+    {
+      LANG: Fortran
+    } "Set psi for divergence cleaning initially to zero"
+  }
+  if (track_divB)
+  {
+    schedule GRHydro_DivBInit IN HydroBase_Initial
+    {
+      LANG: Fortran
+    } "Set divB initially to zero"
+  }
 }
 
 #################################################
diff --git a/src/GRHydro_CalcUpdate.F90 b/src/GRHydro_CalcUpdate.F90
index 81de7db..452b070 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
+  CCTK_REAL :: idx, alp_l, alp_r, Bvec_l, Bvec_r
 
   CCTK_INT :: type_bits, atmosphere, not_atmosphere
   
@@ -55,7 +55,7 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
 
     if (use_weighted_fluxes == 0) then
 
-      !$OMP PARALLEL DO PRIVATE(i,j,itracer,alp_l, alp_r)
+      !$OMP PARALLEL DO PRIVATE(i,j,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
@@ -95,6 +95,13 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
                        (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
                        alp_r * psidcflux(i,j,k)) * idx
                endif
+               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 
+               endif
             endif
     
             if (evolve_tracer .ne. 0) then
@@ -248,6 +255,13 @@ subroutine UpdateCalculation(CCTK_ARGUMENTS)
                      (psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
                      psidcflux(i,j,k)) * idx
              endif
+             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) + ( Bvec_l - Bvec_r ) * idx 
+             endif
           endif
           
         enddo
diff --git a/src/GRHydro_CalcUpdateM.F90 b/src/GRHydro_CalcUpdateM.F90
deleted file mode 100644
index 9bfdb75..0000000
--- a/src/GRHydro_CalcUpdateM.F90
+++ /dev/null
@@ -1,255 +0,0 @@
-  /*@@
-   @file      GRHydro_CalcUpdateM.F90
-   @date      Oct 10, 2010
-   @author    Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
-   @desc 
-   Calculates the update terms given the fluxes. Moved to here so that 
-   @enddesc 
- @@*/
-
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-#include "cctk_Functions.h"
-#include "SpaceMask.h"
-
- /*@@
-   @routine    UpdateCalculationM 
-   @date       Oct 10, 2010
-   @author     Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke
-   @desc 
-   Calculates the update terms from the fluxes.
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   Moved out of the Riemann solver routines to make the FishEye /
-   weighted flux calculation easier.
-   @endhistory 
-
-@@*/
-
-
-subroutine UpdateCalculationM(CCTK_ARGUMENTS)
-    
-  implicit none
-
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_FUNCTIONS
-
-  CCTK_INT :: i,j,k,itracer
-  CCTK_REAL :: idx, alp_l, alp_r
-
-  CCTK_INT :: type_bits, atmosphere, not_atmosphere
-  
-  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
-  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere",&
-                              "in_atmosphere")
-  call SpaceMask_GetStateBits(not_atmosphere, "Hydro_Atmosphere",&
-                              "not_in_atmosphere")
-
-  idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction)
-
-  if (CCTK_EQUALS(method_type, "RSA FV")) then
-
-    if (use_weighted_fluxes == 0) then
-
-      !$OMP PARALLEL DO PRIVATE(i,j,itracer,alp_l, alp_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
-            
-            alp_l = 0.5d0 * (alp(i,j,k) + &
-                 alp(i-xoffset,j-yoffset,k-zoffset))
-            alp_r = 0.5d0 * (alp(i,j,k) + &
-                 alp(i+xoffset,j+yoffset,k+zoffset))
-            
-            densrhs(i,j,k) = densrhs(i,j,k) + &
-                 (alp_l * densflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * densflux(i,j,k)) * idx
-            srhs(i,j,k,1) = srhs(i,j,k,1) + &
-                 (alp_l * sxflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * sxflux(i,j,k)) * idx
-            srhs(i,j,k,2) = srhs(i,j,k,2) + &
-                 (alp_l * syflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * syflux(i,j,k)) * idx
-            srhs(i,j,k,3) = srhs(i,j,k,3) + &
-                 (alp_l * szflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * szflux(i,j,k)) * idx
-            taurhs(i,j,k) = taurhs(i,j,k) + &
-                 (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * tauflux(i,j,k)) * idx
-            Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + &
-                 (alp_l * Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * Bvecxflux(i,j,k)) * idx
-            Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + &
-                 (alp_l * Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * Bvecyflux(i,j,k)) * idx
-            Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + &
-                 (alp_l * Bveczflux(i-xoffset,j-yoffset,k-zoffset) - &
-                 alp_r * Bveczflux(i,j,k)) * idx
-            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 (evolve_tracer .ne. 0) then
-               do itracer=1,number_of_tracers
-                  cons_tracerrhs(i,j,k,itracer) = cons_tracerrhs(i,j,k,itracer) + &
-                       (alp_l * cons_tracerflux(i-xoffset,j-yoffset,k-zoffset,itracer) - &
-                       alp_r * cons_tracerflux(i,j,k,itracer)) * idx
-               enddo
-            end if
-
-            if (evolve_Y_e .ne. 0) then
-               Y_e_con_rhs(i,j,k) = Y_e_con_rhs(i,j,k) + &
-                    (alp_l * Y_e_con_flux(i-xoffset,j-yoffset,k-zoffset) - &
-                    alp_r * Y_e_con_flux(i,j,k)) * idx
-            end if
-            
-            if (wk_atmosphere .eq. 1) then
-
-              if ( (atmosphere_mask(i,j,k) .eq. 1) .or. &
-                   (SpaceMask_CheckStateBitsF90(space_mask,i,j,k,type_bits,atmosphere)) ) then
-
-!!$                We are in the atmosphere so the momentum flux must vanish
-
-                srhs(i,j,k,:) = 0.d0
-
-                if ( ( (atmosphere_mask(i-1,j  ,k  ) .eq. 1) .and. &
-                       (atmosphere_mask(i+1,j  ,k  ) .eq. 1) .and. &
-                       (atmosphere_mask(i  ,j-1,k  ) .eq. 1) .and. &
-                       (atmosphere_mask(i  ,j+1,k  ) .eq. 1) .and. &
-                       (atmosphere_mask(i  ,j  ,k-1) .eq. 1) .and. &
-                       (atmosphere_mask(i  ,j  ,k+1) .eq. 1) &
-                     ) .or. &
-                     ( (SpaceMask_CheckStateBitsF90(space_mask,i-1,j  ,k  ,type_bits,atmosphere)) .and. &
-                       (SpaceMask_CheckStateBitsF90(space_mask,i+1,j  ,k  ,type_bits,atmosphere)) .and. &
-                       (SpaceMask_CheckStateBitsF90(space_mask,i  ,j-1,k  ,type_bits,atmosphere)) .and. &
-                       (SpaceMask_CheckStateBitsF90(space_mask,i  ,j+1,k  ,type_bits,atmosphere)) .and. &
-                       (SpaceMask_CheckStateBitsF90(space_mask,i  ,j  ,k-1,type_bits,atmosphere)) .and. &
-                       (SpaceMask_CheckStateBitsF90(space_mask,i  ,j  ,k+1,type_bits,atmosphere)) & 
-                     ) &
-                    ) then
-
-!!$                    All neighbours are also atmosphere so all rhs vanish
-
-                    densrhs(i,j,k) = 0.d0
-                    taurhs(i,j,k)  = 0.d0
-!!$
-!!$ We should still evolve the B-field in the atmosphere
-!!$
-
-                end if
-              end if
-
-            end if
-
-          enddo
-        enddo
-      enddo
-      !$OMP END PARALLEL DO
-      
-    else
-      
-      call CCTK_WARN(0, "Not supported")
-      
-!!$    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
-!!$          
-!!$          alp_l = 0.5d0 * (alp(i,j,k) + &
-!!$               alp(i-xoffset,j-yoffset,k-zoffset))
-!!$          alp_r = 0.5d0 * (alp(i,j,k) + &
-!!$               alp(i+xoffset,j+yoffset,k+zoffset))
-!!$          
-!!$          densrhs(i,j,k) = densrhs(i,j,k) + &
-!!$               (alp_l * &
-!!$               &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$               &densflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$               alp_r * &
-!!$               &cell_surface(i,j,k,flux_direction) * &
-!!$               &densflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$          sxrhs(i,j,k) = sxrhs(i,j,k) + &
-!!$               (alp_l * &
-!!$               &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$               &sxflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$               alp_r * &
-!!$               &cell_surface(i,j,k,flux_direction) * &
-!!$               &sxflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$          syrhs(i,j,k) = syrhs(i,j,k) + &
-!!$               (alp_l * &
-!!$               &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$               &syflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$               alp_r * &
-!!$               &cell_surface(i,j,k,flux_direction) * &
-!!$               &syflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$          szrhs(i,j,k) = szrhs(i,j,k) + &
-!!$               (alp_l * &
-!!$               &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$               &szflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$               alp_r * &
-!!$               &cell_surface(i,j,k,flux_direction) * &
-!!$               &szflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$          taurhs(i,j,k) = taurhs(i,j,k) + &
-!!$               (alp_l * &
-!!$               &cell_surface(i-xoffset,j-yoffset,k-zoffset,flux_direction) * &
-!!$               &tauflux(i-xoffset,j-yoffset,k-zoffset) - &
-!!$               alp_r * &
-!!$               &cell_surface(i,j,k,flux_direction) * &
-!!$               &tauflux(i,j,k)) * idx / cell_volume(i,j,k)
-!!$          
-!!$        enddo
-!!$      enddo
-!!$    enddo
-
-    end if
-  
-  else if (CCTK_EQUALS(method_type, "Flux split FD")) then
-
-    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
-          
-          densrhs(i,j,k) = densrhs(i,j,k) + &
-               (densflux(i-xoffset,j-yoffset,k-zoffset) - &
-                densflux(i,j,k)) * idx
-          srhs(i,j,k,1) = srhs(i,j,k,1) + &
-               (sxflux(i-xoffset,j-yoffset,k-zoffset) - &
-                sxflux(i,j,k)) * idx
-          srhs(i,j,k,2) = srhs(i,j,k,2) + &
-               (syflux(i-xoffset,j-yoffset,k-zoffset) - &
-                syflux(i,j,k)) * idx
-          srhs(i,j,k,3) = srhs(i,j,k,3) + &
-               (szflux(i-xoffset,j-yoffset,k-zoffset) - &
-                szflux(i,j,k)) * idx
-          taurhs(i,j,k) = taurhs(i,j,k) + &
-               (tauflux(i-xoffset,j-yoffset,k-zoffset) - &
-               tauflux(i,j,k)) * idx
-          Bvecrhs(i,j,k,1) = Bvecrhs(i,j,k,1) + &
-               (Bvecxflux(i-xoffset,j-yoffset,k-zoffset) - &
-               Bvecxflux(i,j,k)) * idx
-          Bvecrhs(i,j,k,2) = Bvecrhs(i,j,k,2) + &
-               (Bvecyflux(i-xoffset,j-yoffset,k-zoffset) - &
-               Bvecyflux(i,j,k)) * idx
-          Bvecrhs(i,j,k,3) = Bvecrhs(i,j,k,3) + &
-               (Bveczflux(i-xoffset,j-yoffset,k-zoffset) - &
-               Bveczflux(i,j,k)) * idx
-          if(clean_divergence.ne.0) then
-             psidcrhs(i,j,k) = psidcrhs(i,j,k) + &
-                  (psidcflux(i-xoffset,j-yoffset,k-zoffset) - &
-                  psidcflux(i,j,k)) * idx
-          endif
-          
-        enddo
-      enddo
-    enddo
-      
-  end if
-  
-  return
-
-end subroutine UpdateCalculationM
-  
diff --git a/src/GRHydro_FluxM.F90 b/src/GRHydro_FluxM.F90
index 77bee7a..fadbfe6 100644
--- a/src/GRHydro_FluxM.F90
+++ b/src/GRHydro_FluxM.F90
@@ -47,21 +47,10 @@ subroutine num_x_fluxM(dens,sx,sy,sz,tau,Bx,By,Bz,&
 !!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w
   tauf = tau*vxt + psipstar*velm - ab0*psiBx/w
 
-!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k]
+!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k)]
   bxf = 0.0
   byf = psiBy * vxt - psiBx*vyt
   bzf = psiBz * vxt - psiBx*vzt
 
 end subroutine num_x_fluxM
 
-
-
-
-
-
-
-
-
-
-
-
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90
index c5a15bb..56e19a1 100644
--- a/src/GRHydro_HLLEM.F90
+++ b/src/GRHydro_HLLEM.F90
@@ -172,6 +172,11 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
 
         avg_det =  SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh)
 
+        call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
+             avg_det,gxxh, gxyh, gxzh, &
+             gyyh, gyzh, gzzh)
+          
+
         vxtp = prim_p(2)-avg_betax/avg_alp
         vytp = prim_p(3)-avg_betay/avg_alp
         vztp = prim_p(4)-avg_betaz/avg_alp
@@ -218,7 +223,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
                  vxtp,vytp,vztp,pressstarp,bxlowp,bylowp,bzlowp,ab0p,wp, &
                  avg_det,avg_alp,avg_beta)
             if(clean_divergence.ne.0) then
-               psidcf = cons_p(6)
+               f1(6)=f1(6)+uxxh*psidcp
+               f1(7)=f1(7)+uxyh*psidcp
+               f1(8)=f1(8)+uxzh*psidcp
+               psidcf = ch_dc*ch_dc*cons_p(6)
             endif
 
           else if (flux_direction == 2) then
@@ -228,7 +236,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
                  vytp,vztp,vxtp,pressstarp,bylowp,bzlowp,bxlowp,ab0p,wp, &
                  avg_det,avg_alp,avg_beta)
             if(clean_divergence.ne.0) then
-               psidcf = cons_p(7)
+               f1(6)=f1(6)+uxyh*psidcp
+               f1(7)=f1(7)+uyyh*psidcp
+               f1(8)=f1(8)+uyzh*psidcp
+               psidcf = ch_dc*ch_dc*cons_p(7)
             endif
 
           else if (flux_direction == 3) then
@@ -238,7 +249,10 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
                  vztp,vxtp,vytp,pressstarp,bzlowp,bxlowp,bylowp,ab0p,wp, &
                  avg_det,avg_alp,avg_beta)
             if(clean_divergence.ne.0) then
-               psidcf = cons_p(8)
+               f1(6)=f1(6)+uxzh*psidcp
+               f1(7)=f1(7)+uyzh*psidcp
+               f1(8)=f1(8)+uzzh*psidcp
+               psidcf = ch_dc*ch_dc*cons_p(8)
             endif
             
           else
@@ -247,10 +261,6 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
           
         else !!! The end of this branch is right at the bottom of the routine
             
-          call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, &
-               avg_det,gxxh, gxyh, gxzh, &
-               gyyh, gyzh, gzzh)
-          
           if (flux_direction == 1) then
             usendh = uxxh
           else if (flux_direction == 2) then
@@ -302,8 +312,14 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
                  avg_det,avg_alp,avg_beta)
 
             if(clean_divergence.ne.0) then
-               psidcfp = cons_p(6)
-               psidcfm = cons_m(6)
+               fminus(6)=fminus(6)+uxxh*psidcm
+               fminus(7)=fminus(7)+uxyh*psidcm
+               fminus(8)=fminus(8)+uxzh*psidcm
+               fplus(6)=fplus(6)+uxxh*psidcp
+               fplus(7)=fplus(7)+uxyh*psidcp
+               fplus(8)=fplus(8)+uxzh*psidcp
+               psidcfp = ch_dc*ch_dc*cons_p(6)
+               psidcfm = ch_dc*ch_dc*cons_m(6)
             endif
 
           else if (flux_direction == 2) then
@@ -329,8 +345,14 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
                  vytm,vztm,vxtm,pressstarm,bylowm,bzlowm,bxlowm,ab0m,wm, &
                  avg_det,avg_alp,avg_beta)
             if(clean_divergence.ne.0) then
-               psidcfp = cons_p(7)
-               psidcfm = cons_m(7)
+               fminus(6)=fminus(6)+uxyh*psidcm
+               fminus(7)=fminus(7)+uyyh*psidcm
+               fminus(8)=fminus(8)+uyzh*psidcm
+               fplus(6)=fplus(6)+uxyh*psidcp
+               fplus(7)=fplus(7)+uyyh*psidcp
+               fplus(8)=fplus(8)+uyzh*psidcp
+               psidcfp = ch_dc*ch_dc*cons_p(7)
+               psidcfm = ch_dc*ch_dc*cons_m(7)
             endif
 
           else if (flux_direction == 3) then
@@ -356,8 +378,14 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
                  vztm,vxtm,vytm,pressstarm,bzlowm,bxlowm,bylowm,ab0m,wm, &
                  avg_det,avg_alp,avg_beta)
             if(clean_divergence.ne.0) then
-               psidcfp = cons_p(8)
-               psidcfm = cons_m(8)
+               fminus(6)=fminus(6)+uxzh*psidcm
+               fminus(7)=fminus(7)+uyzh*psidcm
+               fminus(8)=fminus(8)+uzzh*psidcm
+               fplus(6)=fplus(6)+uxzh*psidcp
+               fplus(7)=fplus(7)+uyzh*psidcp
+               fplus(8)=fplus(8)+uzzh*psidcp
+               psidcfp = ch_dc*ch_dc*cons_p(8)
+               psidcfm = ch_dc*ch_dc*cons_m(8)
             endif
 
           else
@@ -389,8 +417,16 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS)
 
           if(clean_divergence.ne.0) then
              psidcdiff = psidcm - psidcp
-             psidcf = (charmax * psidcfp - charmin * psidcfm + &
-                  charmax * charmin * psidcdiff) / charpm
+
+!!$             psidcf = (charmax * psidcfp - charmin * psidcfm + &
+!!$                  charmax * charmin * psidcdiff) / charpm
+
+!!$   Wavespeeds for psidc are +/-c, not Fast Magnetosonic?
+
+             psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + &
+                  1.d0 * (-1.d0) * psidcdiff) / 2.d0
+
+
           endif
 
         end if !!! The end of the SpaceMask check for a trivial RP.
diff --git a/src/GRHydro_PreLoop.F90 b/src/GRHydro_PreLoop.F90
index bfd4f2b..e0b3bc2 100644
--- a/src/GRHydro_PreLoop.F90
+++ b/src/GRHydro_PreLoop.F90
@@ -94,3 +94,32 @@ subroutine GRHydro_Scalar_Setup(CCTK_ARGUMENTS)
 
 
 end subroutine GRHydro_Scalar_Setup
+
+
+
+ /*@@
+   @routine    GRHydro_DivBInit
+   @date       Apr 06, 2011
+   @author     Bruno Mundim
+   @desc 
+   Set divB=0 initially.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+subroutine GRHydro_DivBInit(CCTK_ARGUMENTS)
+  
+  implicit none
+  
+  DECLARE_CCTK_ARGUMENTS
+
+  divB=0.0
+
+end subroutine GRHydro_DivBInit
+
+
diff --git a/src/GRHydro_SourceM.F90 b/src/GRHydro_SourceM.F90
index 2a5de76..257a24c 100644
--- a/src/GRHydro_SourceM.F90
+++ b/src/GRHydro_SourceM.F90
@@ -116,6 +116,10 @@ subroutine SourceTermsM(CCTK_ARGUMENTS)
   if (clean_divergence .ne. 0) then
      psidcrhs=0.d0
   endif
+
+  if (track_divB .ne. 0) then
+     divB=0.d0
+  endif
   
 
 !!$  Set up the array for checking the order. We switch to second order
@@ -295,22 +299,12 @@ subroutine SourceTermsM(CCTK_ARGUMENTS)
           dy_alp = DIFF_Y_2(alp)
           dz_alp = DIFF_Z_2(alp)
 
-          if(clean_divergence.ne.0) then
-             dx_psidc = DIFF_X_2(psidc)
-             dy_psidc = DIFF_Y_2(psidc)
-             dz_psidc = DIFF_Z_2(psidc)
-          endif
-
         else
 
           dx_alp = DIFF_X_4(alp)
           dy_alp = DIFF_Y_4(alp)
           dz_alp = DIFF_Z_4(alp)
-          if(clean_divergence.ne.0) then
-             dx_psidc = DIFF_X_4(psidc)
-             dy_psidc = DIFF_Y_4(psidc)
-             dz_psidc = DIFF_Z_4(psidc)
-          endif
+
         end if
         
         if (local_spatial_order .eq. 2) then
@@ -436,21 +430,12 @@ subroutine SourceTermsM(CCTK_ARGUMENTS)
              rhohstarW2*invalp*(vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz) -&
              bt*(bxlow*dz_betax + bylow*dz_betay + bzlow*dz_betaz)
 
-        if(clean_divergence.ne.0) then
-           bvcx_source = uxx*dx_psidc+uxy*dy_psidc+uxz*dz_psidc
-           bvcy_source = uxy*dx_psidc+uyy*dy_psidc+uyz*dz_psidc
-           bvcz_source = uxz*dx_psidc+uyz*dy_psidc+uzz*dz_psidc
-        endif
-
         densrhs(i,j,k) = 0.d0
         srhs(i,j,k,1)  = alp(i,j,k)*sqrtdet*sx_source
         srhs(i,j,k,2)  = alp(i,j,k)*sqrtdet*sy_source
         srhs(i,j,k,3)  = alp(i,j,k)*sqrtdet*sz_source
         taurhs(i,j,k) = alp(i,j,k)*sqrtdet*tau_source
         if(clean_divergence.ne.0) then
-           Bvecrhs(i,j,k,1)  = alp(i,j,k)*sqrtdet*bvcx_source
-           Bvecrhs(i,j,k,2)  = alp(i,j,k)*sqrtdet*bvcy_source
-           Bvecrhs(i,j,k,3)  = alp(i,j,k)*sqrtdet*bvcz_source
            psidcrhs(i,j,k) = -1.d0*ch_dc*ch_dc/cp_dc/cp_dc*psidc(i,j,k)
         endif